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Abstract 

Motivated  by  a  factory  scheduling  problem,  we  consider  the  problem  of  input  control 
(subject  to  a  specified  input  mix)  and  priority  sequencing  in  a  multistation,  multiclass 
queueing  network  with  general  service  time  distributions  and  a  general  routing  structure. 
The  objective  is  to  minimize  the  long-run  expected  average  number  of  customers  in  the 
system  subject  to  a  constraint  on  the  long-run  expected  average  output  rate.     Under 
balanced  heavy  loading  conditions,  this  scheduling  problem  can  be  approximated  by  a 
control  problem  involving  Brownian  motion.    Linear  programming  is  used  to  reduce  the 
workload  formulation  of  this  control  problem  to  a  constrained  singular  control  problem 
for  a  multidimensional  Brownian  motion.     The  finite  difference  approximation  method 
is  then  used  to  find  a  linear  programming  solution  to  the  latter  problem.    The  solution 
is  interpreted  in  terms  of  the  original  queueing  system  in  order  to  obtain  an  effective 
scheduHng  policy.  The  priority  sequencing  policy  is  based  on  dynamic  reduced  costs  from 
a  Unear  program,  and  the  workload  regulating  input  policy  releases  a  customer  into  the 
system  whenever  the  workload  process  enters  a  particular  region.  An  example  is  provided 
that  illustrates  the  procedure  and  demonstrates  its  effectiveness. 

Subject  ciassification;  Production/scheduling:  priority  sequencing  in  a  stochastic  job 
shop.  Queues:  Brownian  models  of  network  scheduhng  problems. 
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Harrison  [7]  has  introduced  a  Brownian  system  model  that  approximates  a  multiclass 
queueing  network  with  dynamic  scheduling  capability.  Under  balanced  heavy  loading 
conditions,  this  model  allows  a  queueing  network  scheduling  problem  to  be  approximated 
by  a  control  problem  involving  Brownian  motion.  In  Wein  [20],  a  particular  Brownian 
control  problem  was  solved  that  approximates  the  problem  of  input  control  (subject  to  a 
specified  product  mix)  and  priority  sequencing  in  a  two-station  multiclass  queueing  network 
with  general  service  time  distributions  and  a  general  routing  structure.  The  objective  was 
to  minimize  the  long-run  expected  average  number  of  customers  in  the  system  subject  to 
a  constraint  on  the  long-run  expected  average  output  rate.  The  solution  to  the  Brownian 
control  problem  was  interpreted  in  terms  of  the  original  queueing  system  in  Wein  [21]  in 
order  to  obtain  an  effective  input  control  and  priority  sequencing  policy. 

In  this  paper  we  extend  these  results  from  the  setting  of  a  two-station  network  to  a 
network  with  any  finite  number  of  stations.  The  two-station  Brownian  control  problem 
was  solved  in  [20]  by  (1)  reformulating  the  problem  in  terms  of  workloads,  (2)  using  linear 
programming  to  reduce  the  workload  formulation  to  a  constrained  singular  control  problem 
for  a  one-dimensional  Brownian  motion,  and  (3)  finding  a  closed-form  solution  to  the 
constrained  singular  control  problem.  The  resulting  priority  sequencing  policy  was  based 
on  dynamic  reduced  costs  from  the  linear  program,  and  the  input  policy  depended  on  a 
two-dimensional  workload  process,  which  measured  the  total  expected  amount  of  work  in 
the  network  for  each  of  the  two  stations.  This  worA-ioad  regulating  input  policy  released 
a  job  into  the  network  whenever  the  workload  process  entered  a  particular  region  in  the 


nonnegative  orthant  of  R^;   this  region  was  based  on  the  solutions  to  both  the  hnear 
program  and  the  constrained  singular  control  problem. 

For  the  general  multistation  problem  considered  here,  the  Brownian  control  problem 
can  again  be  reformulated  in  terms  of  workloads  ajid  linear  programming  can  be  employed 
to  reduce  the  workload  formulation  to  a  constrained  singular  control  problem.  Therefore, 
the  resulting  sequencing  policy  is  again  based  on  dynamic  reduced  costs  derived  from  the 
hnear  program.  However,  the  constrained  singular  control  problem  now  involves  a  multi- 
dimensional Brownian  motion  process,  and  the  problem  appezirs  to  be  very  difficult  to 
solve  in  closed  form.  Instead,  we  employ  the  method  of  finite  difference  approximations 
(see  Kushner  [12]  for  a  detailed  development)  to  obtain  a  numerical  solution  to  the  con- 
strained control  problem.  By  discretizing  both  state  and  time,  this  technique  allows  us  to 
approximate  a  controlled  diffusion  (and  functionals  of  the  controlled  diffusion)  by  a  con- 
trolled Markov  chain  (and  functionals  of  the  controlled  Markov  chain).  In  particular,  if  we 
ignore  the  constraints  in  our  constrained  singular  control  problem,  then  the  problem  can 
be  approximated  by  a  controlled  Markov  chain  with  a  long-run  average  cost  criterion.  It 
is  well-known  (see,  for  example,  Manne  [15]  and  Derman  [4])  that  the  latter  problem  can 
be  formulated  and  solved  as  a  linear  program.  Moreover,  we  show  that  under  the  finite 
difference  approximation,  the  constraints  in  our  singular  control  problem  become  linear 
constraints  in  the  linear  programming  formulation  of  the  Markov  chain  control  problem. 
Therefore,  these  constraints  can  simply  be  added  to  the  linear  program  for  the  Markov 
chain  control  problem,  and  an  approximate  solution  to  the  constrained  singular  control 
problem  can  be  found  by  solving  a  linear  program.  As  in  Wein  [21],  the  proposed  input 
policy  is  a  workload  regulating  input  policy,  where  the  region  of  release  depends  on  the 
solution  to  a  linear  program  and  the  approximate  solution  to  the  constrained  singular 
control  problem. 

In  order  to  rigorously  justify  the  finite  difference  approximation,  one  needs  to  prove 
(see  Kushner  [12]-[13])  that  the  optimally  controlled  Markov  chain  (suitably  interpolated) 
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converges  to  the  optimally  controlled  diffusion,  cind  that  the  optimal  cost  of  the  controlled 
Markov  chciin  converges  to  the  optimzJ  cost  of  the  constrained  singular  control  problem. 
Such  a  justification,  which  is  typically  based  on  weak  convergence  methods,  is  not  at- 
tempted here.  However,  we  do  show  that  this  finite  diff"erence  approximation  technique, 
when  applied  to  the  numerical  example  of  a  two-station  queueing  network  scheduling  prob- 
lem considered  in  Wein  [21],  agrees  with  the  closed  form  solution  of  this  problem  derived 
in  Wein  [20]. 

As  in  Wein  [21],  the  customer  release  policy  and  the  priority  sequencing  policy  derived 
here  work  together  in  a  coordinated  way.  The  input  policy  reluctantly  pulls  work  into  the 
network,  in  that  a  customer  is  released  whenever  a  server  is  threatened  with  idleness  and 
there  is  not  too  much  work  already  present  in  the  network.  The  priority  sequencing  policy 
aggressively  pushes  this  work  to  the  various  stations  in  order  to  avoid  server  idleness. 
As  will  be  seen  in  Section  10,  the  proposed  policies  offer  a  significant  improvement  in 
performance  over  conventional  customer  release  and  priority  sequencing  policies. 

This  paper  is  organized  as  follows.  In  Section  1,  the  queueing  network  scheduling 
problem  is  described,  and  the  workload  formulation  of  the  approximating  Brownian  con- 
trol problem  is  given  in  Section  2.  Linear  programming  is  used  in  Section  3  to  reduce  the 
workload  formulation  to  a  constrained  singular  control  problem,  which  is  described  in  Sec- 
tion 4.  The  finite  difference  approximation,  which  allows  the  constrained  singular  control 
problem  to  be  approximated  by  a  constrained  Markov  chain  control  problem,  is  given  in 
Section  5,  and  a  linear  programming  solution  to  the  latter  problem  is  derived  in  Section  6. 
In  Sections  7  and  8,  respectively,  the  solution  to  the  workload  formulation  is  interpreted  in 
terms  of  the  original  queueing  system  in  order  to  obtain  effective  priority  sequencing  and 
input  control  policies,  respectively,  for  the  original  queueing  network  scheduling  problem. 
In  Section  9,  a  heuristic  extension  to  the  input  policy  is  described  that  allows  for  the  de- 
cision of  which  class  of  customer  to  next  release  into  the  system,  not  just  when  to  release 
the  customer.  An  example  is  given  in  Section  10  that  illustrates  the  entire  procedirre  and 
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demonstrates  its  effectiveness,  aiid  concluding  remarks  are  offered  in  Section  11. 

1.     The  Queueing  Network  Scheduling  Problem 

The  queueing  network  scheduling  problem  studied  in  this  paper  is  motivated  by  a 
scheduling  problem  that  is  encountered  in  many  factories;  see  Wein  [21]  for  a  description 
of  the  motivating  factory  scheduling  problem.  Consider  a  queueing  network  with  I  single- 
server  stations  and  A'  customer  classes.  Class  k  customers  require  service  at  a  pcirticular 
station  s{k)  and  have  service  times  that  are  independent  and  identically  distributed  random 
variables  with  mean  m  jt  and  variance  s]..  A  class  k  customer,  upon  completion  of  service  at 
station  s(k),  turns  next  into  a  class  j  customer  with  probability  Pkj  and  exits  the  network 
with  probability  1  —  Ylj=i  Pkj-  The  Markovian  switching  matrix  P  has  spectral  radius  less 
than  one,  and  so  all  customers  eventually  exit  the  network  with  probability  one.  Because 
the  number  of  classes  is  allowed  to  be  arbitrary,  this  routing  structure  is  almost  perfectly 
general. 

There  are  input  control  and  priority  sequencing  decisions  in  the  scheduling  problem. 
We  assume  there  is  an  ample  supply  of  customers  who  are  waiting  to  gain  entry  into  the 
network,  and  that  each  of  these  customers  has  an  exogenously  specified  class  designation. 
The  class  designations  are  such  that  qk  is  the  long-run  proportion  of  class  k  customers 
released  into  the  network.  The  entering  class  mix  vector  q  =  {qk)  is  such  that  ^;(.=i  9^  ~ 
1.  We  will  assume  that  the  class  designations  are  deterministic;  however,  they  could  be 
Markovian  without  changing  the  nature  of  the  analysis.  The  input  decisions  are  to  choose 
the  non-decreasing  process  A'^  =  {N{t),t  >  0},  where  N(t)  is  the  cumulative  number  of 
customers  released  into  the  network  in  the  time  inteval  [0,<].  Thus,  the  input  decisions 
allow  for  full  discretion  over  the  timing  of  the  release  of  customers  into  the  system,  but 
do  not  allow  for  the  choice  of  which  class  of  customer  to  release.  In  Section  9,  we  develop 


a  heuristic  scheme  that  allows  the  controller  to  decide  which  class  of  customer  to  release 
into  the  system;  this  scheme  guarantees  that  the  actual  mix  that  is  released  is  sufficiently 
close  to  the  desired  mix  q. 

The  sequencing  decisions  are  to  dynamically  choose  which  class  of  customer  to  serve 
at  each  station  in  the  network.  Although  preemptive  resume  scheduling  is  assumed,  the 
assumptions  regarding  preemption  do  not  have  an  impact  on  the  proposed  scheduling 
policy. 

There  is  a  per  unit  time  holding  cost  Cfc  that  is  incurred  when  a  class  k  customer  is 
in  the  queueing  network.  Define  the  throughput  rate  of  the  queueing  network  to  be  the 
number  of  customer  departures  per  unit  of  time.  Then  our  queueing  network  scheduling 
problem  is  to  find  the  input  and  sequencing  policy  that  minimizes  the  long-run  expected 
average  holding  costs  subject  to  a  constraint  that  the  long-run  expected  average  throughput 
rate  is  greater  than  or  equal  to  the  specified  lower  bound  A.  If  ct  =  c  for  k  =  l,...,/\, 
then  the  objective  is  to  minimize  the  long-run  expected  average  number  of  customers  in 
the  network.  Since  the  throughput  constraint  will  be  tight  in  general,  this  latter  objective 
is  equivalent  (by  Little's  formula)  to  minimizing  the  long-run  expected  average  cycle  time 
of  customers,  where  a  customer's  cycle  time  is  the  length  of  time  it  spends  in  the  network. 

2.     The  Workload  Formulation  of  the  Brownian  Control  Problem 

Harrison  [7]  has  shown  how  the  queueing  network  scheduling  problem  described  in  the 
last  section  is  approximated  by  a  control  problem  for  a  Brownian  network.  In  Wein  [20],  it 
is  shown  how  this  Brownian  control  problem  is  reformulated  in  terms  of  workloads.  Since 
the  proposed  scheduling  policy  depends  only  on  the  solution  to  the  workload  formulation, 
we  will  go  directly  to  the  workload  formulation  of  the  Brownian  control  problem  that 
approximates  the  queueing  network  scheduHng  problem  described  in  Section  1. 
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Let  p  =  {pi)  be  the  /—vector  of  server  utilizations,  or  traffic  intensities,  for  the  /  sta- 
tions; later  in  this  section,  p  will  be  explicitly  defined  in  terms  of  the  problem  parameters. 
The  Brownian  approximation  assumes  the  existence  of  a  large  integer  n  such  that 

\/n{\  —  Pi)  is  positive  and  of  moderate  size  for  i  =  1, ...,  /.  (2-1) 

Let  Qk  =  {Qk{t)J  ^  0}  be  the  number  of  class  k  customers  in  the  system  at  time  t,  for 
k  =  1,...,A',  and  let  /,  =  {I,{t),t  >  0}  be  the  cumulative  amount  of  time  that  the  server 
at  station  i  is  idle  in  the  time  interval  [0,<],  for  i  =  1, ...,/. 

With  the  parajTieter  n  fixed,  define  the  scaled  queue  length  process  Zk  =  {Zk{t),t  >  0} 

by 

Z^m  =  ^4^,  /  >  0  and  ^^  =  1,...,  a;  (2.2) 

and  the  scaled  cumuiative  idleness  process  U,  =  {Ut{t),t  >  0}  by 

U,{t)  =  ^^^,  t>0  and  1=^1,. ..J.  (2.3) 

Define  the  one-dimensional  scaled  centered  input  process  6  by 

Xnt-  N(nt) 
m  = ^^-^,  i>0,  (2.4) 

where  N{t)  is  the  cumulative  number  of  customers  released  into  the  system  up  to  time  t 
and  A  is  the  specified  average  throughput  rate.  The  processes  Z  =  {Zk),U  —  (U,),  and  0 
axe  the  control  processes  in  the  workload  formulation  of  the  Brownian  control  problem. 
Let  the  AT— vector  A  =  (A^t)  de  defined  by 

A  =  q~X.  (2.5) 

Since  q  is  the  entering  class  mix  vector,  it  follows  that  A^  represents  the  average  number 
of  class  k  customers  that  must  depart  from  the  system  per  unit  of  time  in  order  to  satisfy 
the  throughput  rate  constraint.  Define  the  A'  x  A'  input-output  matrix  R  =  (Rkj)  by 

Bkj=m-\6,k-Pjk\  (2.6) 
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where  8jk  =  \  \i  j  —  k  and  6jjt  =  0  otherwise.  The  term  Rkj  represents  the  average  rate 
at  which  class  k  customers  axe  depleted  when  class  j  customers  are  being  served.  Since 
the  routing  matrix  P  is  transient,  the  matrix  R  is  nonsingular  and  there  exists  a  unique 
nonnegative  solution  /?  =  (/3t)  to  the  flow  balance  equations 

A  =  i?/?.  (2.7) 

Define  The  I  x  K  resource  consumption  matrix  A  =  (Anc)  by 

\  0,     otherwise. 
Then  the  server  utilization  vector  p  referred  to  in  condition  (2.1)  is  defined  by 

p  =  AI3.  (2.9) 

Now  define  the  /  x  A'  worA-ioad  profile  matrix  M  =  (Mjt)  by 

M  =  AR-\  (2.10) 

The  element  Mik  represents  the  total  expected  remaining  amount  of  work  for  a  class 
k  customer  at  station  i  until  the  customer  exits  the  network.  Let  the  /—dimensional 
worA'ioad  process  W  be  defined  by 

W{t)  =  MZ{t),  t  >  0,  (2.11) 

so  that  Wi{t)  is  interpreted  as  the  total  expected  amount  of  scaled  work  anywhere  in  the 
network  for  station  i  at  time  t.  Now  define  the  /—vector  v  =  {v,)  by 

V  =  Mq,  (2.12) 

so  that  Vi  is  interpreted  as  the  expected  total  amount  of  time  over  the  long-run  that  server 
i  spends  on  each  customer.  By  Proposition  1  in  Wein  [21], 

/9,  =i),Afor  i  =  1,...,/.  (2.13) 
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Let  the  /—vector  7  =  (7,)  be  defined  by 

-,,  =   y/^{l-p,)i0Tl  =   l,...,I.  (2.14) 

Let  C{i)  be  the  set  of  all  customer  classes  k  such  that  s{k)  =  i,  and  define  the  i^— vector 
a  =  (Qk)  by 

ak  =  —  for  alike  C{i).  (2.15) 

Now  define  X  to  be  a  A'— dimensional  Brownian  motion  process  with  drift  6  and  covairiance 
E,  where 

6  =  y/Ti^X  -  Ra)  (2.16) 

and 

K 

^ji  =  I][c^^"^ir'^itj(<5j7  -  Pki)  +  Qkm-'slRjkRik].  (2.17) 

Finally,  let  B  be  the  /—dimensional  Brownian  motion  process  defined  by 

B{t)  =  MXit),  t  >  0.  (2.18) 

The  process  B  has  drift  M6  and  covariance  M'EM'^.  It  was  shown  in  Wein  [21]  that  the 
drift  vector  MS  =  — y,  where  7  was  defined  in  (2.14). 

The  Brownian  control  problem  is  obtained  by  letting  the  system  parameter  n  defined 
in  the  balanced  heavy  loading  condition  (2.1)  approach  infinity.  Exactly  the  same  notation 
used  for  the  scaled  processes  Z,  U,  and  6  are  used  in  defining  the  Brownian  control  problem 
in  order  to  retain  the  queueing  network  interpretation  of  the  Brownian  model.  Define  the 
workload  formulation  of  the  Brownian  control  problem  as  choosing  right  continuous  with 
left  limit  (RCLL)  processes  Z,  U  and  6  (K-,  I-  and  one-dimensional,  respectively)  so  eis  to 

1  /^   ^ 

minimize    limsup— £r[/      >     CkZk{t)dt]  (2.19) 

subject  to    Z,  U    and  6  are  non  —  anticipating  with  respect  to  A',  (2.20) 
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U  is  non  —  decreasing  with  U{0)  =  0,  (2-21) 

Z{t)  >   0   for  all  t  >  0,  (2.22) 

limsup-£;[l/,(r)]    <    7.    for  i  =  1, ...,/,    and  (2.23) 

T— oo     T 

MZ{t)  =  B{t)  +  U{t)  -  ve{t)   for  all   t  >  0.  (2.24) 

Problem  (2.19)-(2.24)  is  referred  to  as  the  workload  formulation  because  the  basic 
system  state  equation  (2.24)  is  in  terms  of  the  /—dimensional  workload  process  W  defined 
in  (2.11). 


3.     Solving  For  Z  in  Terms  of  U 

In  this  section,  we  express  the  optimal  process  Z  in  the  workload  formulation  in  terms 
of  the  control  process  U .  Suppose  we  are  given  a  process  U  that  satisfies  constraints  (2.20), 
(2.21),  and  (2.23).  Then  the  optimal  Z  and  0  processes  are  found  by  solving  the  following 
linear  program  at  each  time  t: 

K 

min       Y  CkZk{i)  (3.1) 

K 
subject  to    ^  A/a-^Jt(<)  +  v,e{t)  =  B,it)  +  Ihit)  for  i  =  1,...,/,  (3.2) 


fc=i 


Zkit)  >   0,    for  k  =  1,...,A'.  (3.3) 


At  each  time  t  >  0,  this  linear  program  may  have  a  different  set  of  right  hand  side 
values.  The  dual  of  this  linear  program  will  be  easier  to  analyze,  since  it  has  a  static 
constraint  set.  Let  us  define  the  dual  variables  7r(i)  =  (7r,(i))  and  state  the  duai  linear 
program  to  be  solved  at  time  t: 

I 

max         y[B,{t)  +  U,{t)]n,{t)  (3.4) 

1=1 
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subject  to    y_]  M,>7ri(<)  <    c^    for  k  —  \,  ...,K,  (3.5) 

1=1 
/ 
J2vini{t)  =  0.  (3.6) 

1=1 

Before  analyzing  the  dual  LP  (3.4)-(3.6),  let  us  define  the  (/— l)-dimensional  workload 
imbalajice  process  W  =  {Wt{t))  by 

Wi(t)  =  piW,{t)  -  p,Wiit),  t  >  0,  for  i  =  1, ...,  7-1.  (3.7) 

By  (2.11),  (2.13),  (2.24)  and  (3.7),  the  workload  imbalance  process  can  also  be  expressed 
as 

W,(t)  =  pjB,(t)  -  p,Bi(i)  +  pjU,{i)  -  p,Uiitl    t  >  0,  for  i  =  1, ...,  7-1.  (3.8) 

Thus,  the  workload  imbalance  process  does  not  depend  on  the  control  process  6,  and  is 
expressed  in  terms  of  the  Brownian  motion  process  B  and  the  control  process  U. 

Returning  to  the  dual  LP  (3.4)-(3.6),  use  (3.6)  to  eliminate  7r/(t)  from  the  problem  amd 
use  (2.13)  to  substitute  the  utilization  levels  p  for  the  vector  v  to  obtain  the  (7  —  1)- variable 
dual  linear  program 

/-I 
max         p7^^W^.(07r.(0  (3.9) 

jr,(0,...,T/_i(()  ^ 

:=1 

/-I 

subject  to   c^^  y^{piM,k  —  PtMjk)T^i{t)  <   pj    for  k  =  1,...,  A'.  (3.10) 

1=1 

Denote  the  solution  to  (3.9)-(3.10)  by  (7r*(<), ...,  7rJ_j(<))  and  solve  for  7rJ(f)  using 
equation  (3.6).  Let  Ck{t)  denote  the  dynamic  reduced  cost  for  the  variable  Zk{t)  in  the 
linear  program  (3.1)-(3.3).  That  is, 

c;t(0  =  c, -^7r'(f)A/.t.  (3.11) 

1=1 
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We  will  return  to  the  dynamic  reduced  costs  in  Section  4,  where  the  costs  become  the 
basis  for  our  proposed  sequencing  policy.  By  the  analysis  above,  it  follows  that  the  linear 
program  (3.1)- (3. 3)  can  be  expressed  as 

min    Y^CkZkit)  (3.12) 

K 

subject  to    ^{piMik  -  piMik)Zk{t)  =  Wi{t)  for  i  =  1, ...,  /  -  1,  (3.13) 

it=i 

Zk{t)>   0,    for  Jfc  =  l,...,A".  (3.14) 

Denoting  the  solution  to  this  linear  program  by  Z^(t),  we  construct  the  optimal  queue 
length  process  Z*  from  Zl{i),k  =  1,...,A',  for  all  /  >  0.  Notice  that  this  optimal  process 
does  not  depend  on  the  control  process  0  and  depends  on  the  control  process  U  only 
through  the  workload  imbalance  process  W  defined  in  (3.7). 

4.     The  Resulting  Control  Problem 

In  this  section  we  substitute  the  optimal  queue  length  process  Z*  for  Z  into  the 

workload  formulation  (2.19)-(2.24),  and  reduce  this  problem  to  one  of  finding  the  optimal 

/—dimensional  cumulative  idleness  process  U .    By  duality  theory,  it  is  known  that  the 

optimal  value  of  the  primal  and  dual  objectives  in  problems  (3.1)-(3.3)  and  (3.9)-(3.10) 

will  be  equal,  or  that 

K  I-\ 

Y^ckZi{i)  =  p-;'Y.w,{t)7^:{t).  (4.1) 

Define  the  function  h:  R''^  -^  R^  by 


k=\  «=1 


7-1 

h{W{t))  =  pj''£w,it)n:it).  (4.2) 

t=i 

Then  h  is  a  piecewise-linear,  continuous  function  of  W{t),  and  h  achieves  a  minimum  of 
zero  at  the  point  W{t)  —  0. 
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Define  the  (/  —  1)— dimensional  Brownian  motion  B  by 


B,it)  =  p,B,(t)  -  p,Bj(i),  <  >  0,  for  i  =  1, ...,  /  -  1. 


(4.3) 


Recalling  that  the  two-dimensioneJ  Brownian  motion  B  has  drift  —7,  it  follows  that  B  has 
drift  ft  =  {fi,),  where 

Mi  =  V^ipi  -  Pl),  for  t  =  1, ...,/-  1.  (4.4) 

Notice  that  when  the  system  is  perfectly  balanced  (i.e.,  p,  =  p  for  i  =  1,...,/),  then  the 
drift  /z  =  (0,  ...,0).  Define  the  (/  -  1)  x  /  matrix  T  by 


T  = 


/PI      0      0 
0^/0 


0       0      0 
V  0       0      0 


0 


-PI    \ 

-p2 


PI        0        -pi-2 

0     pi     -pi-\  / 


(4.5) 


Then  B  has  (/  -  1)  X  (/  -  1)  covariance  matrix  a  =  TMT,M'^T'^. 

The  resulting  control  problem  is  to  find  a  non-anticipating  (with  respect  to  the 
A'— dimensional  Brownian  motion  A'),  non-decreasing,  /—dimensional,  RCLL  process  U 
to 


minimize    lim  sup  —  Ej- 

T—'oo     ■'■ 


I    h{W{t))dt 
Jo 


subject  to    hmsup —Ei[Ui(T)]    <   7,,  for  i  =  1, ...,/, 


W,{t)  =  B,{t)  +  piUi{t)  -  p,Ui{t),  for  I  =  1, ...,  /  -  1. 


(4.6) 

(4.7) 
(4.8) 


In  problem  (4.6)-(4.8),  the  controller  observes  the  (/  —  1)— dimensional  Brownian 
motion  process  B  and  exerts  the  non-anticipating  /—dimensional  control  U  =  (Ui).  The 
constraint  (4.7)  determines  an  upper  bound  on  the  long-run  expected  average  amount  of 
control  exerted.  The  resulting  controlled  process  is  the  (/  —  1)  — dimensional  workload 
imbalance  process  W,  and  since  the  optimal  process  Z*  derived  in  the  previous  section  is 
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used,  costs  axe  incurred  at  a  rate  h(W{t)).  Notice  that  for  i  =  1, ...,  /  —  1,  the  control  Ui 
only  affects  the  ith  component  of  W,  whereas  Uj  affects  all  /  —  1  components  of  W. 

Problem  (4.6)-(4.8)  will  be  referred  to  as  a  constrained  singular  control  problem. 
The  name  "singular"  stems  from  the  fact  that  the  state  of  the  controlled  process  can 
be  instantaneously  changed  by  the  controller  and,  as  a  result,  the  optimal  process  U  is 
continuous  but  singular  (that  is,  the  set  of  time  points  at  which  U  increases  has  measure 
zero).  Problem  (4.6)-(4.8)  can  be  viewed  as  a  variant  of  the  Unite-fuel  problem  for  Brownian 
motion,  in  which  there  is  a  constraint  on  the  toted  cimount  of  controlling  effort  that  can  be 
exerted  (or  fuel  that  can  be  consumed).  In  the  traditional  finite-fuel  problem,  the  amount 
of  control  exerted  is  constrained  to  be  finite  over  some  finite  or  infinite  time  interval;  see 
Benes,  Shepp,  and  Witsenhausen  [1],  Chow,  Menaldi,  and  Robin  [3],  and  Kaxatzas  [10]  for 
variants  of  this  problem.  In  contrast,  the  constraints  (4.7)  are  on  the  long-run  expected 
average  amount  of  control  exerted. 

If  a  solution  U*  to  the  constrained  control  problem  (4.6)-(4.8)  can  be  found,  then  the 

optimal  9  process  is,  via  equation  (3.2), 

K 

e'{i)  =  v;'[B,{t)  +  U;ii)  -Y,MuZ:{t)]   for  aU  i  >  0,  (4.9) 

where  Z*{t)  is  the  solution  to  (3.12)-(3.14).  The  optimal  process  6*  is  never  explicitly  used 
to  develop  the  proposed  sequencing  and  input  control  policies. 

We  now  show  that  the  cost  function  h  in  the  constrained  singular  control  problem  is 
convex. 

Proposition  4.1.    The  function  h{W{t))  defined  in  (4.2)  is  convex  in  W{t). 

Proof.    By  definition, 

7-1 

h{w)  =      max      pj^  /^  Wi^^i  (4-10) 

7ri,...,7r;_i  ^—f 


t=l 


subject  to 


7-1 


<^k^^ipi^^tk-piMjk)n,  <   PI   for  k  =  1,...,A'.  (4.11) 


1=1 
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Let  tt"  be  the  solution  to 

7-1 


max      p7^Vu"f7r,  (4.12) 

ri,...,jr;_i  ^-— ' 


1=1 
,6 


subject  to  (4.11),  and  let  ir°  be  the  solution  to 

7-1 

mzix     pj^y^  w^iTi  (413) 

jri,...,ir/_i  ^-^ 

1=1 

subject  to  (4.11).  Thus,  hiw")  =  pj^  ^.1=1  "">?  and  h{w'')  =  p/'  ^fj/  u'^Trf.  For  some 
A  €  [0, 1],  let  w^  =  (1  —  A)u'°  +  Au'*,  and  suppose  n^  is  the  solution  to 

7-1 

max      pJ^y^w^TTi  (4.14) 

Tl,...,7r7_i  •t-—' 

1=1 

subject  to  (4.11);  thus,  h{w^)  =  pj^  E[=/[(1  "  ^)^?  +  AtifJTr,^.  Therefore 

7-1  7-1 

(1  -  X)h{w'')  +  Xhiiv'')  =  (1  -  A)p7^  J]u.°<  +  Ap7^  ^u-fTrf  (4.15) 

1=1  1=1 

7-1  7-1 

>  (1  -  A)p7^  J]u.>,^  +  Ap7^  J^zifTT.^  (4.16) 

1=1  i=l 

=  Mti'),  (4.17) 


where  inequality  (4.16)  holds  because  n^  satisfies  constraint  (4.11).    | 

The  goal  of  the  next  two  sections  is  to  find  a  solution  to  the  constrained  singular 
control  problem.  To  that  end,  we  will  make  one  more  minor  transformation  of  problem 
(4.6)-(4.8).  The  goal  of  this  transformation  is  to  ehminate  the  coefficients  in  front  of  the 
Ui  terms  in  equation  (4.8).  In  the  perfectly  balanced  case  where  p,  =  p  for  i  =  1,...,/, 
the  workload  imbalance  process  W  and  the  Brownian  motion  B  could  have  been  defined 
by  W,{t)  =  W,{t)  -  Wj{t)  =  E*=,(^^.;t  -  Mik)ZUt)  and  B,(t)  =  B,{t)  -  B;(0  for 
i  =  1, ...,  /  —  1;  then  equation  (4.8)  could  be  expressed  as  Wi{t)  =  B,{t)  +  U,{t)  —  Ui{t). 
When  the  system  is  not  perfectly  balanced,  let  pi  be  defined  by 

7-1 
Px=     n     />;forz  =  l,...,/-l.  (4.18) 
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and  malce  the  following  definitions: 


I  I 

U:{t)=     n     PjUi{t)  Bind  ^:  =     J]     p,7.  fori  =  !,...,/; 


(4.19) 


W°{t)  =  piW,{t)  and  B-{t)  =  p,B,(t)  for  i  =  1, ...,  /  -  1. 


(4.20) 


Also,  define  the  function  h°  by  h°{W°(t))  =  h{W(t)).  Notice  that  h"  will  be  a  piecewise 
linear,  convex,  and  continuous  function  with  a  minimum  of  zero  at  zero.  Then  problem 
(4.6)-(4.8)  can  be  expressed  as  choosing  a  non-decreasing,  non-anticipative  (with  respect 
to  X),  /—dimensional,  RCLL  process  U°  to 

■r 


...       ,.  1  ^ 

mimmize    limsup— r/j- 
T— oo    T 


1 


/    h\W\t))dt 
Jo 


subject  to    limsup —Ei[U°{T)]    <   7°,  for  i  =  1, ...,  J, 
T->oo    T 


W°{t)  =  BUt)  +  U°{t)  -  Unt),  for  I  =  1, ...,  7-1. 


(4.21) 

(4.22) 
(4.23) 


For  ease  of  notation,  we  shall  drop  all  of  the  "0"  superscripts  and  subsequently  analyze 
the  following  problem:  choose  a  non-decreasing,  non-anticipating  (with  respect  to  X), 
/—dimensional,  RCLL  process  U  to 


1  f  '. 

minimize    limsup— £'r     /     h[\\  {t))dt 


subject  to    \\msup —Ej:[Ui{T)]    <   7,,  for  i  =  1,...,/, 


W,(t)  =  B,{t)  +  U,{t)  -  Ui{t),  for  i  =  1,...,/-  1. 


(4.24) 

(4.25) 
(4.26) 


The  drift  and  covariance  of  the  Brownian  motion  B  will  be  denoted  by  the  (/—  1)— dimen- 
sional vector  //  =  (fj.,)  and  the  (7  —  1)  x  (/  —  1)  matrix  a  =  (a,j),  respectively. 

5.     The  Finite  Difference  Approximation  Method 
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One  possible  approach  to  analyzing  the  constrained  singular  control  problem  is  to 
form  the  Lagrangian  relaxation  of  (4.24)-(4.26),  where  the  constraints  (4.25)  axe  placed  in 
the  objective  function  and  multiplied  by  the  (unknown)  Lagrange  multipliers  /  =  (/,): 

/ 


min  limsup  —  iTr 


I  h{w{t))dt  +  yi,u,{T) 


(5.1) 


subject  to   W,{t)  =  B,{t)  +  U,{t)  -  Ui{t),  for  i  =  1, ...,  7-1.  (5.2) 

Under  the  case  where  the  multipliers  /  are  known,  the  Lagrangian  problem  has  been  am- 
aJyzed  by  Taksar  [19]  for  the  case  /  =  2  and  by  Meneddi  et  al.  [16]  for  general  /.  The 
optimality  conditions  for  the  Lagrangian  problem  ase  to  find  {V,g,U)  such  that 

i^Tr  ^T/  .^T/ 

Min  {TV(x)  +  h{x)-gj,  +  -—,...,/;_,  + , /,  _  Y^         }  =  o,  (5.3) 

where 

/-I  1-1  ^2  /-I         ^ 

1=1  J=l  '      J       ,=1  • 

is  the  generator  of  B.  In  (5.3),  the  potential  function  V(x)  :  R'~^  — >  R^  is  the  cost 
incurred  under  the  optimal  policy  when  the  initial  state  of  the  controlled  process  ly  is  x 
minus  the  cost  incurred  under  the  optimal  policy  when  the  initial  state  of  W  is  zero.  Also, 
the  gain  g  appearing  in  (5.3)  is  the  minimal  average  cost  per  unit  time,  independent  of 
initial  state. 

An  argument  identical  to  Theorem  5.1  in  Wein  [20]  shows  that  the  constrained  problem 
(4.24)-(4.26)  can  be  solved  by  finding  a  solution  (including  the  approximate  multipliers  /) 
to  (5.3)-(5.4)  that  simultaneously  satisfies  constraints  (4.25).  For  the  special  case  /  =  2, 
a  closed  form  solution  (V,g,  /j ,  l2,Ui ,  U2)  to  (5.3)-(5.4)  and  (4.25)  was  found  in  Wein  [20]. 
Although  Menaldi  et  al.  [IG]  have  shown  that  (5.3)-(5.4)  are  sufficient  conditions  (along 
with  some  regularity  conditions)  for  optimality  of  the  Lagrangian  problem  (5.1)-(5.2)  and 
have  proven  the  existence  of  a  solution  {V,g,U)  to  (5.3)-(5.4),  there  appears  to  be  little 
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hope  of  finding  a  closed  form  solution  to  (5.3)-(5.4).  Therefore,  finding  a  closed  form 
solution  {V,g,l,U)  to  (5.3)-(5.4)  and  (4.25)  does  not  seem  possible. 

Thus,  our  goal  is  to  obtain  a  numerical  solution  to  the  constrained  singular  control 
problem.  However,  finding  a  numerical  solution  {V,g,  /,  U)  also  appears  to  be  an  arduous 
task,  since,  even  if  a  numerical  solution  {V,g,  U)  to  (5.3)-(5.4)  were  obtained  for  some  set 
of  multipliers  /  (not  an  obvious  task  in  itself),  a  search  over  the  space  of  /  would  be  needed 
to  assure  that  U  satisfied  constraints  (4.25).  However,  for  the  case  1  =  2,  Wein  [20]  showed 
that  the  solution  (V,gJi,l2,U\,  U2)  included  any  nonnegative  pair  (/i ,  /a)  such  that  li  plus 
I2  equaled  a  particular  constant.  If  a  solution  {V,g,  I,  U)  to  (5.3)-(5.4),  (4.25)  included  any 
nonnegative  multipliers  (/i, ...,  //)  such  that  ^,_j  U  equaled  a  particular  constant  (I  do  not 
know  if  this  is  true),  then  a  reduction  in  the  space  of  /  would  be  possible. 

Considering  these  difficulties,  a  better  approach  to  a  numerical  solution  to  this  prob- 
lem appears  to  be  the  finite  difference  approximation  technique  developed  by  Kushner  [12]. 
The  basic  idea  behind  this  method  is  to  systematically  approximate  a  controlled  diffusion 
process  by  a  controlled  finite  state  Markov  chain,  and  then  to  solve  the  resulting  optimiza- 
tion problem  for  the  controlled  Markov  chain.  Weak  convergence  methods  are  then  used 
to  verify  that  the  controlled  diffusion  process  (and  its  corresponding  optimal  cost)  can  be 
approximated  arbitrarily  closely  by  the  controlled  Markov  chain  (and  its  optimal  cost).  In 
this  paper,  the  finite  difference  method  is  described  and  two  numerical  examples  are  given, 
but  no  weaic  convergence  proof  is  provided. 

This  method  is  particularly  powerful  for  the  control  problem  (4.24)-(4.26)  because  of 
the  existence  of  constraints  (4.25).  Although  very  little  is  known  about  multi-dimensional 
stochastic  control  problems  with  side  constraints  using  a  dynamic  programming  formula- 
tion, the  finite  difference  approximation  can  easily  incorporate  the  side  constraints  (4.25), 
as  will  be  shown  in  the  next  section. 

We  now  begin  to  describe  the  finite  difference  method  and,  for  ease  of  reference,  we 
will  retain  most  of  the  notation  of  Kushner  [12].  When  there  are  /  stations  in  the  network, 
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the  state  space  of  W  in  the  control  problem  is  R^~^ .  In  order  to  numerically  solve  the 
problem,  we  need  to  confine  the  process  W  to  a  bounded  set,  which  will  be  denoted  by  G. 
The  set  G  needs  to  be  large  enough  so  that  W  never  reaches  the  boundary  of  G,  which 
is  denoted  by  dG.  Define  h,  not  to  be  confused  with  the  cost  function  in  (4.24),  to  be 
the  finite  difference  interval,  which  dictates  how  finely  the  state  space  R^~^  is  discretized. 
Let  i?[~'  be  the  finite  difference  grid  on  R^~^\  a  point  x  €  -R/,"^  if  there  exists  integers 
ni,...,n/_i  such  that  x  =  X^,Ji  heiUi,  where  e,  is  the  unit  vector  in  the  tth  coordinate 
direction.  The  approximating  controlled  Markov  chain,  which  we  denote  by  {^^,n  >  0}, 
will  have  state  space  G/,  =  Rf,~^  H  G.  Before  we  define  the  transition  probabilities  for  the 
controlled  Markov  chain,  the  controls,  or  actions,  will  be  described. 

Let  U{x)  be  the  action  set  for  the  controlled  Markov  chain  {^^,n  >  0}  when  it  is  in 
state  X  G  Gh-  Recall  that  the  actions  in  the  Markov  chain  control  problem  correspond  to 
the  controls  U  in  the  singular  control  problem.  Let  u^  =  1  if  the  control  corresponding  to 
Ui  is  exerted,  and  let  u^  =  0  if  the  control  corresponding  to  Ui  is  not  exerted,  for  i  =  !,...,/. 
Then  action  u  G  U{x)  is  defined  by  u  =  (uj,  ...,u/),  and  the  cardinality  of  the  action  set 
is  2  for  £l11  X  ^  dG.  We  will  not  concern  ourselves  with  the  action  set  and  transition 
probabilities  when  x  €  dG,  since,  in  practice,  the  controlled  Markov  chain  will  never  reach 
the  boundary;  one  just  enlcirges  the  set  G  if  the  controlled  Markov  chain  reaches  dG. 

Let  P  {x,y;u)  denote  the  transition  probability  from  state  x  to  state  y  when  the 
action  u  =  (ui,...,u/)  is  used  in  state  x  G  Gh-  Define 

/-I  7-1  7-1 

and  assume  that 

7-1 

a„  -     Y^     |a,j|  >  Ofor  z  =  1,...,/- 1.  (5.6) 

As  mentioned  on  page  92  of  Kushner  [12],  if  condition  (5.G)  does  not  hold,  then  a  trans- 
formation to  the  principal  vectors  of  a  can  be  applied  to  assure  (5.6)  in  a  new  coordinate 
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system.  By  definition  of  u,  the  action  u  —  (0,...,0)  corresponds  to  exerting  no  control.  In 
this  case,  define  the  transition  probabilities  P^[x,y;u)  by 


>/.,„    .^„    ^...,_«..-E0:;#.}l«.;l+2/^/^. 


P\x,x±  e,h;u)  =  '"':,'         -^  (5-'^) 

af. 
P  {x,x  +  Cih  +  e,/?;  u)  =  P(x,x  —  Cih  —  tjh;  u)  =      '"'    for  i  ^  j,  (5-8) 

a~ 

P  (x,  x  —  Cih  +  ejh)  =  P{x,  x  +  Cih  —  tjh)  =  — ^  for  i  ^  j,  (5-9) 

and  P'^{x,y;u)  =  0  otherwise.     Notice  that  the  transition  probabilities  P''{x,y;u)  are 

nonnegative  and  sum  over  y  to  one  for  each  x. 

When  the  action  u  ^  (0,...,0),  then  the  Markov  chain  transitions  are  deterministic, 

due  to  the  exertion  of  controls.  In  this  case, 

/-I  7-1 

P^{x,  x  +  Y^  e,u,h  -  Y^  eiuih;  u)  =  1  (5.10) 

1=1  :=1 

and  P'^{x,y;u)  =  0  otherwise. 

Define  the  interpolation  interval  At    by 

Af"  =  -^.  (5.11) 

Equation  (5.11)  relates  the  size  of  the  discretized  time  intervals  to  the  size  of  the  discretized 
space  intervals  and,  together  with  the  transition  probabihties  P''{x,y]u)  in  (5.7)-(5.10), 
assures  that  the  first-  and  second-order  moments  of  i^^i  —  (^,  conditioned  on  ^^,  are 
consistent  with  those  of  the  controlled  diffusion  process  W. 

Thus,  the  controlled  process  W  described  by  equation  (4.26)  has  been  approximated 
by  the  controlled  Markov  chain  {^^,n  >  0}  with  transition  probabilities  given  by  (5.7)- 
(5.10).  Since  there  are  no  costs  on  the  controls  exerted  in  problem  (4.24)-(4.26),  the  cost 
incurred  when  action  u  is  taken  and  the  controlled  Markov  chain  is  in  state  x  is  simply  given 
by  h{x),  where  the  function  h,  which  appears  in  objective  (4.24),  is  defined  in  (4.2).  Thus, 
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ignoring  constraints  (4.25)  for  the  moment,  our  problem  (4.24)-(4.26)  is  approximated 
by  a  problem  of  controlling  a  finite-state,  finite  action-set  Markov  chain  with  a  long-run 
average  cost  criterion.  In  the  next  section,  we  show  how  to  express  the  constraints  (4.25)  in 
terms  of  the  approximating  Markov  chain  control  problem,  ajid  how  to  solve  the  resulting 
constrained  Markov  chEiin  control  problem. 

6.     A  Linear  Programming  Solution 

It  is  well  known  that  a  Markov  chain  control  problem  with  a  long-run  average  cost 
criterion  can  be  formulated  and  solved  as  a  linear  program;  readers  are  referred  to  Manne 
[15]  and  Dermaji  [4]  for  early  work  on  this  subject.  We  begin  this  section  by  formulating 
the  controlled  Markov  chain  problem  described  in  the  last  section,  which  ignored  the 
constraints  (4.25),  as  a  linear  program.  Suppose  that  the  cardinality  of  the  set  G/,  is  .A/, 
and  let  the  states  be  indexed  by  x  =  l,...,Af  and  the  actions  indexed  by  u  =  1,...,2^.  A 
stationary  policy  for  the  Markov  chain  control  problem  will  be  described  by  /?  =  {0x{u),  x  = 
1, ...,  A/;  u  =  1,  ...,2^),  where  0z{u)  equals  the  probability  that  action  u  is  chosen  when  the 
controlled  Markov  chain  is  in  state  x.  If  ^i(u)  equals  zero  or  one  for  all  x  =  1, ...,  M  and 
u  =  1,...,2^,  then  ^  is  referred  to  as  a  pure  stationary  policy;  otherwise,  /9  is  referred  to 
as  a  randomized  stationary  policy.  Let  tt  =  (tt^-u)  denote  the  steady-state  probability  that 
the  controlled  Markov  chain  will  be  in  state  x  and  action  u  will  be  chosen  if  policy  /3  is 
used;  thus,  tt  is  policy- dependent  and  must  satisfy 

7r^„  >  Ofor  X  =  1,...,A/;  u  =  l,...,2^  (6.1) 

M      2' 

^  J]7r,„  =  l,  and  (6.2) 

r=l  u=l 
2'  M      2' 

E^y«  =  ^^7r,„P''(x,y;u)for  j/  =  l,...,A/,  (6.3) 

U=l  1  =  1   u  =  l 

where  the  transition  probabilities  P'^{x,y\u)  were  given  in  equations  (5.7)-(5.10).    Since 
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the  expected  average  cost  under  policy  /3  is  5Zj._i  X^„=i  T^zuh{x),  the  Markov  chain  control 

problem  is  to  find  tt  =  (t^zu)  to 

\f    2' 
min  V] /_]  7rj.u/i(x)  (6-4) 

1=1 u=l 

subject  to  constraints  (6.1)-(6.3).  If  tt*  =  (7r*„)  solves  the  linear  program  (6.1)-(6.4),  then 

the  optimal  policy  jS*  is 

^x(")  =  ^J^"       ■  (6.5) 

Thus,  the  linear  program  (6.1)-(6.4)  solves  the  problem  (4.24)-(4.26)  that  ignores 
constraints  (4.25).  We  now  show  that  under  the  finite  difference  approximation,  the  / 
constraints  (4.25)  can  be  expressed  as  /  linear  constraints  on  tt  =  (tTj-u).  From  equation 
(5.10),  if  the  control  u  in  the  controlled  Markov  chain  {^n,  u  >  0}  is  such  that  //,  =  1,  then 
the  corresponding  cumulative  amount  of  control  exerted  by  J7,  in  problem  (4.24)-(4.26) 
is  increased  by  the  finite  difference  interval  size  h.  Equation  (5.10)  also  implies  that  the 
steady-state  probability  that  u,  =  1  is  ^^.-1  X](uu  =1)  ^'«-  Thus,  by  the  finite  difference 
approximation,  the  term  Wm supj^^^  T~^ Ex[Ui{T)]  in  (4.25)  is  approximated  by 

lim   "^^-^^i-=^^"^"  for  .  =  1, ...,  /,  (6.6) 

"-00  nAt  '  ^      ^ 

which  equals 

h 

by  the  critical  relationship  (5.11)  relating  the  time  and  space  intervals.     Therefore  the 

constraints  (4.25)  are  approximated  by 

A/        2'  , 

E      E     ^x«  <77^fort  =  l,...,/.  (6.8) 

Thus,  a  solution  to  the  constrained  Markov  chain  control  problem  and  an  approximate 
solution  to  the  constrained  singular  control  problem  (4.24)-(4.26)  can  be  obtained  by  solv- 
ing the  Hnear  program  (6.1)-(6.4),  (6.8).  The  solution  to  this  linear  program  will  yield  an 
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optimal  policy  /3  for  which  the  controls  for  up  to  /  states  are  rcindomized.  Since  the  cost 
function  h  is  convex  and  achieves  a  minimum  at  zero,  the  optimal  policy  will  be  charac- 
terized by  a  bounded  region  containing  the  origin.  This  region  is  such  that  u  =  (0,  ...,0) 
inside  this  region  and  u  ^  (0,  ...,0)  outside  of  this  region.  The  boundcir}'  of  the  region  acts 
as  a  reflecting  barrier  that  keeps  the  controlled  process  within  the  region  containing  the 
origin.  The  bounded  region  will  play  a  key  role  in  defining  the  customer  release  policy  in 
Section  8. 

The  linear  program  (6.1)-(6.4),  (6.8)  suffers  from  the  curse  of  dimensionality  that  is 
so  often  encountered  in  control  problems.  For  example,  if  the  grid  Gh  is  such  that  each 
coordinate  dimension  is  discretized  into  L  —  l  segments,  then  the  cardinality  of  Gh  is  L^~^ 
and  the  linear  program  has  2^ L^~^  variables  and  L^~^  +/  +  1  constraints,  not  including  the 
nonnegativity  constraints  (6.1).  The  size  of  this  problem  can  be  reduced  in  several  ways. 
Since  the  holding  costs  are  convex  and  have  a  minimum  at  zero,  any  control  u  -^  (0,  ...,0) 
that  pushes  the  controlled  process  further  from  the  origin  will  clearly  be  suboptimal  and 
need  not  be  explicitly  considered  in  the  linear  program.  Also,  the  linear  program  can  be 
solved  several  times  (using  smaller  values  of  the  finite  difference  interval  h  for  later  runs), 
and  portions  of  the  state  space  that  the  controlled  process  never  reaches  can  be  eliminated 
from  the  subsequent  linear  programs,  as  long  as  the  controlled  process  never  reaches  the 
boundcury  of  G/,.  Moreover,  the  finite  difference  interval  h  can  be  state-dependent,  so  that 
a  finer  grid  can  be  used  in  the  sensitive  portion  of  the  grid  Gh  (where  the  optimal  control 
is  uncertain)  and  a  courser  grid  cein  be  used  in  the  insensitive  portion  (where  the  optimal 
control  is  known).  Finally,  some  decomposition  approaches  used  in  linear  programming 
may  be  used  to  exploit  the  structure  of  the  constraint  set  (in  particular,  the  structure  in 
(5.7)-(5.10)  that  appears  in  (6.3));  readers  are  referred  to  Kushner  and  Chen  [14]  for  work 
on  this  topic. 

We  finish  this  section  with  a  numerical  example  that  illustrates  the  accuracy  of  the 
finite  difference  approximation  for  solving  the  constrained  singular  control  problem  (4.24)- 
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(4.26).  The  example  is  taken  from  Section  7  of  Wein  [21]  and  is  derived  from  a  particular 
two-station  queueing  network  scheduling  problem;  the  scheduling  problem  is  identical  to 
the  problem  considered  in  this  paper  with  7  =  2  stations.  The  cost  function  h{x)  in  (4.24) 
is  given  by 

where  hi  =  .101  and  /12  =  .3703.  The  one-dimensional  Brownian  motion  B  has  drift  //  =  0 
and  vEiriance  a  =  10.93.  The  righthand  side  values  in  (4.25)  are  71  =  72  =  -9.  The  closed 
form  solution  to  this  problem  (see  Wein  [20])  is  characterized  by  the  interval  endpoints  I 
and  r,  where 


and 


/  =  -- — —, ^  =  -4.771  (6.10) 

hi  +  h2  271  ^         ' 

r  =  —^;^  =  1.301.  (6.11) 

"1  +  "2  271 

The  controlled  process  W  behaves  as  a  one-dimensional  reflected,  or  regulated,  Brownian 
motion  (see  Harrison  [6]  for  a  detailed  development)  on  the  interval  [—4.771,1.301]  and 
the  optimal  control  processes  U*  and  U2  are  the  local  times  at  the  points  -4.771  and  1.301, 
respectively.  The  optimal  objective  function  value  is  given  by  (see  Wein  [20]) 

t"^\  .  =  .2409.  (6.12) 

47i(/!i+/l2)  ^  ' 

To  formulate  the  linear  program  (6.1)-(6.4),  (6.8),  a  finite  difference  interval  of  size 
h  =  0.05  was  used  and  the  bounded  set  G  was  taken  to  be  the  interval  [—5.0,5.0].  Thus 
the  state  space  is  {—5.0, —4.95, —4.90,  ...,4.95,5.0},  and  the  states  are  indexed  by  x  = 
1,...,201.  The  action  set  U{x)  consists  of  the  four  actions  (0,0),  (1,  0),  (0, 1),  and  (1,1), 
which  correspond,  respectively,  to  exerting  no  control,  exerting  only  Ui,  exerting  only  U2, 
and  exerting  both  Ui  and  U2.  The  non-zero  transition  probabilities  are 

P{x,x  +  1;0,0)  =  P{x,x-  1;0,0)  =  ^,  (6.13) 

P(x,x  +  1;1,0)  =  P{x,x  -  1;0,1}  =  1,  and  (6.14) 

P(x,x;l,l)  =  l.  (6.15) 
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Of  course,  control  u  —  (1, 1)  will  never  be  used  and  control  u  =  (l,0)(u  =  (0, 1),  respec- 
tively) will  never  be  used  when  the  controlled  process  is  positive  (negative,  respectively). 
Since  Qh  =  a  hy  (5.5),  constraints  (6.8)  are  given  by 

^      A  (.05)(.9) 

E      E     ^ru  <^-j^^  .00412  iovr  =  h2.  (6.16) 

1=1  {u:u,  =  l} 

Rather  thcin  indexing  the  states  by  x  =  1,...,201,  let  us  denote  the  state  of  the 
controlled  Markov  chain  by  y  €  [—5.0, —4.95,  ...,4.95,5.0].  The  solution  to  the  linear 
program  yields,  via  equation  (6.5), 

/?;(0,0)  =  1  for  yG  [-4.7,1.25],  (6.17) 

/9l4.8(l,0)  =  l,  (6.18) 

/3:.75(1,0)  =  .327,  (6.19) 

/?;.^5(0,0)  =  .673,and  (6.20) 

^r.3o(0,l)  =  l.  (6.21) 

The  optimal  objective  function  value  for  the  linear  program  was  .2411,  which  is  very  close 
to  the  derived  value  of  .2409  appearing  in  equation  (6.12).  If  we  interpolate  the  policy 
(3  at  y  =  —4.75,  the  approximate  solution  to  the  constrained  singular  control  problem 
is  characterized  by  the  interval  [—4.784,1.300],  where  the  controlled  process  W  behaves 
as  a  one-dimensional  reflected  Brownian  motion  on  this  interval  and  U\  and  U2  are  the 
local  times  of  W  at  the  two  respective  endpoints.  The  interval  derived  from  the  finite 
difference  approximation  [—4.784,1.300]  is  very  close  to  the  exact  interval  [—4.771,1.301] 
at  the  rather  modest  finite  difference  interval  size  of  h  =  .05.  This  particular  example  was 
also  solved  at  the  finer  interval  size  of  /i  =  .01.  The  optimal  objective  function  value  for 
the  linear  program  was  again  .2411,  and  the  interpolated  interval  was  [-4.777,1.300].  Thus, 
the  finite  difference  approximation  method  is  very  accurate,  at  least  for  the  simple  case 
where  a  known  solution  exists. 
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7.  The  Sequencing  Policy 

In  this  section  we  will  interpret  the  solution  to  the  workload  formulation  in  order  to 
obtain  a  priority  sequencing  policy  to  the  original  queueing  network  scheduling  problem. 
As  in  Harrison  [7],  Harrison  and  Wein  [8],  and  Wein  [21],  the  policy  will  be  interpreted  in 
terms  of  the  optimal  control  process  Z*,  where  Zl(t)  is  interpreted  as  the  (scaled)  number 
of  class  k  customers  in  the  system  at  time  t;  readers  are  referred  to  these  previous  works 
for  a  more  detailed  discussion  on  the  interpretation  of  the  solutions  to  Brownian  control 
problems.  In  particular,  the  sequencing  policy  is  based  on  the  dynamic  reduced  costs  Ck{i) 
computed  in  (3.11).  The  reduced  cost  for  a  class  k  customer  at  time  t  is  the  increase  in  the 
optimal  objective  function  value  of  the  linear  program  (3.1)-(3.3)  per  unit  increase  in  the 
righthand  side  of  the  nonnegativity  constraint  Zk{t)  >  0.  This  vcdue  can  be  interpreted  as 
the  extra  cost  incurred  if  one  were  forced  to  hold  a  class  k  customer  in  queue  at  time  t. 
Thus,  the  higher  the  value  of  Ckit\  the  more  expensive  it  is  to  hold  a  class  k  customer  in 
the  system  at  time  t. 

However,  each  customer  class  requires  a  different  amount  of  expected  processing  time 
before  exiting.  Therefore,  it  is  natural  to  consider  the  amount  of  effort  needed  to  completely 
process  a  customer,  in  addition  to  the  cost  incurred  in  holding  the  customer.  As  pointed 
out  by  Yang  [24],  the  ratio 

-^i^lL-  (7.1) 

measures  how  costly  a  class  k  customer  is  at  time  t  per  unit  of  remaining  processing  time. 
Our  proposed  policy  is  to  give  priority  at  each  station  to  the  customer  class  with  the  largest 
value  of  this  dynamic  index.  Yajig  [24]  has  shown  that  this  policy,  when  applied  to  the 
Brownian  analysis  of  a  single-server  multiclass  queue  with  per  unit  time  holding  cost  Ck 
for  class  k  customers  (see  also  Section  6  of  Harrison  [7]),  reduces  to  the  well-known  Cfj.  rule 
(see,  for  example,  Klimov  [11]). 
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Complementan'  slackness  implies  that  the  reduced  cost  Ck{t)  equals  zero  when  Z^{t)  > 
0  in  (3.1)-(3.3).  Thus,  the  policy  proposed  in  (7.1)  will  serve  a  customer  from  classes  with 
Zl{t)  >  0  only  when  there  are  no  customers  present  from  classes  with  Zl(t)  =  0.  Our 
policy  is  therefore  consistent  with  observations  from  existing  heavy  traffic  limit  theorems 
(see,  for  example,  Whitt  [22],  Harrison  [5],  Reiman  [IS],  Johnson  [9],  Peterson  [17],  and 
Chen  eind  Mandelbaum  [2])  that  the  normalized  queue  length  process  of  the  lowest  priority 
customers  is  positive  and  the  normalized  queue  length  processes  of  the  higher  priority 
customers  vanish;  of  course,  these  results  assume  static  (i.e.,  independent  of  the  state 
of  the  queueing  system)  priority  rankings,  whereas  we  are  proposing  dynamic  (i.e,  state- 
dependent)  priority  rankings. 

Notice  that  it  is  possible  for  several  different  customer  classes  with  Zl{t)  >  0  (ajid 
therefore  with  Ck{t)  =  0)  to  be  served  at  a  common  station.  At  times  when  only  these 
customers  are  present  at  a  particular  station,  a  tie-breaking  rule  is  needed  to  decide  which 
of  these  classes  to  serve  next.  We  employ  the  tie-breaking  rule  proposed  in  Yang  [24],  which 
attempts  to  reduce  the  total  expected  holding  cost  X)fc=i  ^kZk{t)  along  a  steepest  descent 
direction.  Pleaders  are  referred  to  that  paper  for  a  derivation  of  this  heuristic  tie-breaking 
rule;  only  a  description  of  the  rule  will  be  presented  here. 

Define  the  (7  —  1)  x  A'  worA'ioad  imbaiance  proiUe  matrix  M  =  {Mik)  by 

M,k  =  piM,k  -  PtMik   for  I  =  1, ...,  I-l;k  =  l, ....  A'.  (7.2) 

Then  constrciint  (3.13)  can  be  expressed  as  J2k=\  ^ikZk{t)  =  W^i(0  fo^  ^  =  li---i-^  ""  1- 
Let  Msit)  be  the  submatrix  of  the  matrix  M  consisting  of  the  columns  that  axe  in  the 
optimal  basis  of  the  linear  program  (3.12)-(3.14)  at  time  t.  It  follows  that  Mfl(/)  is  an 
(/  —  1)  X  (/  —  1)  invertible  matrix  for  all  t.  Let  the  (/  —  1)— dimensional  vector  AZ(/)  be 
defined  by 

AZ{t)  =  MB\t)n'{t),  (7.3) 

where  n*{t)  =  {n^{t),...,n'j_^(t))  is  the  optimal  solution  to  the  dual  LP  (3.9)-(3.10).    If 
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class  k  is  in  the  optimal  basis  at  time  t,  then  we  denote  its  component  of  AZ(t)  in  (7.3) 
by  AZk{t)  in  the  obvious  way.  The  dynamic  tie-breaking  rule  is  to  give  higher  priority  at 
time  t  to  the  classes  with  the  higher  values  of  AZk{t). 

To  summarize  the  priority  sequencing  policy,  we  compute  a  dynamic  reduced  cost 
cic{t)  for  each  class  k  at  each  time  t  using  (3.11).  If  class  k  is  in  the  optimal  basis  at  time 
t,  then  Ck{t)  —  0  and  a  secondary  index  AZjt(i)  is  computed  according  to  equation  (7.3). 
Recall  that  C{i)  is  the  set  of  customer  classes  that  can  be  served  at  station  i.  At  time 
t,  server  i  gives  priority  to  the  customer  class  k  E  C{i)  with  the  largest  value  of  Ck(t).  If 
Ck{t)  =  0  for  all  customers  present  at  station  i  at  time  t,  then  server  i  gives  priority  to  the 
customer  class  with  the  largest  value  of  AZk{t). 

8.  The  Input  Policy 

In  this  section  the  solution  {Z* ,U* ,6*)  to  the  workload  formulation  (2.19)-(2.24)  is 
interpreted  in  order  to  obtain  an  input,  or  customer  release,  policy  for  the  original  queueing 
system.  The  input  policy  will  be  interpreted  in  terms  of  all  three  controls,  in  contrast  to 
the  sequencing  policy,  which  was  interpreted  solely  in  terms  of  the  optimal  queue  length 
process  Z*. 

We  begin  by  interpreting  the  Z*  control.  Notice  that  only  7—1  constraints  in  the 
dual  LP  (3.9)-(3.10)  will  be  tight  at  any  time  <,  and  thus,  by  complementary  slackness, 
only  7—1  components  of  the  optimal  control  process  Z*  can  be  positive  at  any  time  t. 
Therefore,  since  W{t)  =  MZ*{t)  for  all  <  >  0,  the  7— dimensional  workload  process  W 
stays  on  the  boundary  of  a  polyhedral  cone  in  the  nonnegative  orthant  of  R^ .  This  cone, 
which  is  not  necessarily  convex,  is  generated  by  a  number  of  extreme  rays  emanating  from 
the  origin,  where  each  ray  corresponds  to  a  customer  class  that  is  in  the  optimal  basis  of 
the  LP  (3.12)-(3.14)  for  some  value  of  W{t).  Thus  the  number  of  different  extreme  rays, 
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and  hence  the  number  of  different  (/—  1)— dimensional  faces  of  the  cone,  equals  the  number 
of  extreme  points  of  the  dual  constraint  set  (3.10).  See  Figure  1  for  an  example  with  three 
stations  (that  is,  7  =  3)  emd  six  extreme  rays. 


W,  (t) 


FIGURE  1.  W{t)  Stays  on  the  Cone  Boundary. 
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The  workload  process  W  can  also  be  expressed  in  terms  of  the  other  two  controls, 
U*  and  6*.  By  equations  (2.11)  and  (2.24),  the  controller  in  the  workload  formulation 
observes  an  /—dimensional  Brownian  motion  process  B,  exerts  the  optimal  controls  U*, 
which  is  the  /—dimensional  scaled  cumulative  idleness  process,  and  6*,  which  is  the  one- 
dimensional  scaled  centered  input  process,  and  obtains  the  controlled  process  W,  which  is 
the  scaled  workload  process,  via  the  equations 

W,{t)  =  B,{t)  +  U*{t)  -  v,e*{t),  for  ?  =  1, ...,  /  and  f  >  0.  (8.1) 

Recall  that  the  solution  U*  to  the  constrained  singular  control  problem  is  characterized 
by  a  bounded  region  in  R^~^  containing  the  origin.  The  boundary  of  the  region  acts  as  a 
reflecting  barrier  that  keeps  the  (/—  1)— dimensional  workload  imbalance  process  W  within 
this  region  containing  the  origin.  Moreover,  the  control  process  J7*  is  exerted  only  when 
the  workload  imbalance  process  W  reaches  the  reflecting  barrier.  Exerting  the  control  U* 
is  interpreted  as  incurring  server  idleness  at  station  i,  for  i  =  1, ...,/. 

The  restriction  of  the  (/—I )— dimensional  workload  imbalance  process  W"  to  a  bounded 
region  containing  the  origin  implies  the  restriction  of  the  /—dimensional  workload  process 
W  to  a  finite  portion  of  the  cone  boundary.  Thus  the  reflecting  boundary  in  /?  ,  which 
was  derived  in  the  solution  to  the  constrained  singular  control  problem,  essentially  trun- 
cates the  polyhedral  cone  in  R^ .  The  intersection  in  R^  of  the  botmdaxy  of  the  original 
polyhedral  cone  and  the  reflecting  barrier  will  be  referred  to  as  the  upper  edge  of  the 
truncated  cone.  See  Figures  2  and  3  for  typical  cases  when  /  =  2  and  3,  respectively.  In 
Figure  2,  the  reflecting  boimdary  in  R^  is  the  interval  [a,  6],  and  the  upper  edge  of  the 
truncated  cone  consists  of  two  points.  The  upper  edge  of  the  truncated  cone  is  a  curve  in 
R^  in  Figure  3  and,  in  general,  is  of  dimension  1  —  2. 

We  can  now  make  the  following  three  observations  about  the  optimal  solution  to  the 
workload  formulation:  (1)  the  /—dimensional  workload  process  W  stays  on  the  truncated 
cone  boundary,  (2)  the  control  process  U*  is  exerted  only  when  the  workload  process  W 
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W(t) 
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FIGURE  2.  The  Truncated  Cone  in  a  Two-Station  Example. 

reaches  the  upper  edge  of  the  truncated  cone,  and  (3)  when  the  workload  process  is  not  at 
the  upper  edge  of  the  truncated  cone,  then  only  the  input  process  6*  is  used  to  keep  the 
workload  process  on  the  truncated  cone  boundary. 

The  goal  of  this  section  is  to  develop  a  customer  release  policy  for  the  original  queueing 
system  that  operationalizes  these  three  observations.  AU  these  observations  are  developed 
under  the  ideahzed  assumptions  of  the  balanced  heavy  loading  condition  (2.1).  Notice 
that  in  the  actual  queueing  system,  the  workload  process  may  not  reside  exactly  on  the 
truncated  cone  boundary.   The  actual  workload  process  may  reside  in  the  interior  of  the 


30 


W    (t) 
2 


Upper  Edge 


W3(t) 


W.(t) 


FIGURE  3.  The  Truncated  Cone  in  a  Three-Station  Example. 

truncated  cone  or,  since  the  extremal  rays  of  {PF|VF  =  MZ,  Z  >  0}  may  not  coincide  with 
the  extremal  rays  generating  the  truncated  cone,  may  reside  outside  the  truncated  cone. 
We  now  attempt  to  operationalize  the  control  process  9*.  In  order  to  see  how  the  scaled 
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centered  input  process  6*  is  manipulated,  recall  that  by  equation  (2.13)  and  the  balanced 
heavy  loading  condition  (2.1),  v,  is  approximately  equal  to  Vj  for  i,j  =  1,...,/,  and  so  9* 
caji  move  jdong  a  direction  that  is  close  to  the  (1,1,...,1)  direction  in  R^ .  By  definition 
(2.4),  when  6*  is  increased  and  moves  in  the  direction  that  is  close  to  (1, 1, ...,  1),  input  is 
being  increased  relative  to  the  nominal  input  rate  A.  Similarly,  when  6*  is  decreased  emd 
moves  in  the  direction  that  is  close  to  (  — 1,— 1,...,— 1),  input  is  being  witliheld  relative  to 
the  nominal  input  rate. 

When  the  workload  process  W  is  in  the  cone  interior  and  not  near  the  upper  edge, 
then  6*  is  decreased  in  order  to  keep  W  on  the  truncated  cone  boundary.  Similarly,  when 
the  workload  process  W  is  outside  the  cone  and  not  near  the  upper  edge,  then  9*  is 
increased  in  order  to  keep  W  on  the  truncated  cone  boundary.  As  in  Wein  [21],  let  us 
interpret  the  action  "increase  ^*"  to  simply  mean  "release  a  customer  into  the  system" 
cind  the  action  "decrease  ^*"  to  mean  "cease  input".  The  naive  policy  that  emerges  from 
this  interpretation  and  from  observations  (l)-(3)  above  is  to  only  release  a  customer  into 
the  system  at  times  t  when  the  workload  process  W[t)  is  on  the  outside  of  the  truncated 
cone. 

A  precise  definition  of  what  is  meant  by  "outside  the  truncated  cone"  will  be  given 
shortly,  but  first  the  insights  behind  observations  (l)-(3)  will  be  summarized.  When  the 
workload  process  leaves  the  interior  of  the  truncated  cone,  then  the  workload  process 
becomes  too  imbaJanced  (recall  the  importance  of  the  workload  imbalance  process  W)  and 
at  least  one  server  becomes  threatened  with  idleness.  In  this  case,  a  customer  is  released 
into  the  network  to  avoid  server  idleness.  However,  when  the  workload  process  reaches 
the  upper  edge  of  the  truncated  cone,  then  there  is  already  ample  work  in  the  system  and 
the  controller  refuses  to  release  any  more  customers  into  the  system  and  is  willing  to  incur 
server  idleness. 

We  now  return  to  the  issue  of  defining  the  term  "outside  the  truncated  cone".  Let 
Rb  denote  the  bounded  region  in  R^~^  that  characterizes  the  solution  to  the  constrained 
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singular  control  problem  in  Section  6.  Recall  that  the  number  of  faces  on  the  boundary  of 
the  polyhedral  cone  in  R^  is  equal  to  the  number  of  extreme  points  in  the  dual  constraint 
set  (3.10).  Let  us  denote  this  number  by  F,  and  index  the  faces,  and  hence  the  dual 
extreme  points,  by  /  =  1,  ...,F.  For  each  face  /,  let  us  define  Rf  C  R^~^  to  be  the  region 
in  the  workload  imbalance  space  where  extreme  point  /  is  the  optimal  solution  to  the  dual 
LP  (3.9)-(3.10).  For  dual  extreme  point  /,  let  the  customer  classes  in  the  optimal  basis  of 
the  primal  LP  be  indexed  by  A;i, ...,  kj-i.  Then  the  hyperplane  generated  by  cone  face  / 
is  given  by 

Y,aiW,{t)  =  0,  (8.2) 


1=1 
where 


W,{t)      ...      Wi{t) 
Mik,       ...      Mjk, 


(8.3) 


{aiW,{t),...,aiWj{t))  =  dei 

Let  r f  denote  the  halfspace  in  R^  (lying  outside  of  the  cone)  generated  by  the  hyperplane 
(8.2)-(8.3). 

The  workload  regulating  input  policy  will  be  to  release  a  customer  into  the  network 
whenever  the  /—dimensional  workload  process  W(t)  enters  the  release  region  Ui</<F.Ry, 
where 

Rf  =  {W{t)\(W{t)  eRff)  Rb)  n  iW{t)  e  r;  n  R+)},  (8.4) 

and  W{t)  is  defined  in  terms  of  W{t)  in  (3.7);  this  region  defines  the  vague  term  "outside 
the  truncated  cone". 

In  order  to  motivate  policy  (8.4),  it  is  easiest  to  consider  a  two-station  example,  as 
in  Figure  4,  where  the  upper  edge  is  given  by  two  points;  let  us  denote  the  lower  point 
by  (a^i,yi)  and  the  upper  point  by  (^2,2/2)5  where  xi  >  t/i  and  X2  <  y2-  For  the  sake 
of  concreteness  and  without  loss  of  generality,  suppose  the  queueing  network  generating 
Figure  4  was  a  balanced  system,  so  that  Vi  =  V2  and  8*  moves  in  the  (1,1)  direction 
or  the  (—1,-1)  direction;  then  the  refiecting  boundary  [a,b]  introduced  in  Figure  2  is 
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[a'2  —  y2,^i  —  yi]-  Notice  that  as  long  as  the  workload  process  W  is  in  the  darkly  shaded 
region  in  Figure  4  (the  intersection  of  the  darkly  shaded  regions  and  the  gray  regions  are 
along  the  45  degree  Unes),  then  0*  alone  can  be  increased  to  move  W  to  the  truncated 
cone  boundary.  When  W{t)  =  (xi,y),  where  y  <  j/i,  then  U2  alone  is  used  to  move  the 
workload  process  to  the  truncated  cone  bovmdary,  in  this  case  to  the  point  (a;i,yi).  In 
the  lower  (respectively,  upper)  gray  region  of  Figure  2,  some  combination  of  6*  and  U2 
(respectively,  U*)  is  used  to  move  W  to  the  truncated  cone  boundary.  Therefore,  the 
release  region  for  this  example  should  be  at  least  the  darkly  shaded  regions  and  at  most 
the  union  of  the  darkly  shaded  regions  and  the  gray  regions.  Although  either  extreme  (or 
some  compromise  between  the  two  extremes)  would  probably  lead  to  an  effective  customer 
release  policy,  we  have  proposed  in  (8.4)  the  minimvun  release  region,  which  in  this  example 
corresponds  to  the  deirkly  shaded  regions.  This  choice  does  not  maintain  consistency  with 
previous  work  (the  policy  proposed  in  Wein  [21]  was  the  maximal  release  region,  which 
corresponds  to  the  union  of  the  darkly  shaded  regions  and  the  gray  regions  in  Figure  4); 
however,  the  minimal  release  region  suggested  here  is  more  easily  generedizable  to  higher 
dimensions. 

For  the  example  in  Figure  4,  the  reflecting  boundarj'  Rb  is  given  by  W{t)  €  [o,i],  or 
y2  —  X2  <  Wi(0  ~  ^"^2(0  ^  yi  ~  ^1-  If  the  upper  face  (respectively,  lower  face)  is  denoted 
by  face  1  (respectively,  face  2),  then  (see  Wein  [20])  i?i  is  W{t)  <  0,  or  -fj^  <  £i  =  1, 
and  i?2  is  ^.'  ^  >  1.  Recalling  the  definition  of  rf,  it  is  easily  verified  that  the  region 
specified  in  (8.4)  corresponds  to  the  darkly  shaded  regions  in  Figure  4. 

In  order  to  develop  an  explicit  description  of  the  release  region,  an  explicit  expression 
is  required  for  the  (/  —  1)— dimensional  bounded  region  that  characterizes  the  solution 
to  the  constrained  singular  control  problem.  However,  only  a  numerical  solution  to  the 
constrained  singular  control  problem  has  been  derived  in  this  paper.  Thus,  the  numerical 
solution  needs  to  be  transformed  into  an  explicit  expression  for  the  bounded  region.  In 
Section  10,  an  example  is  carried  out  for  a  three-station  network,  and  the  reflecting  bound- 
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FIGURE  4.  The  Release  Region  in  a  Two-Station  Example. 


axy  containing  the  bounded  region  in  R^  is  approximated  by  a  piecewise  hnear  boundary 
in  order  to  develop  an  explicit  release  region  in  R^ .  Presumably,  an  approximation  in  this 
spirit  is  required  to  develop  release  regions  when  /  >  3. 

Notice  that  the  workload  regulating  input  policy  defined  above  would  ignore  a  differ- 
ence that  exists  between  the  actual  queueing  system  and  the  idealized  heavy  traffic  limit. 
As  pointed  out  in  Wein  [21],  this  difference  can  be  understood  by  making  the  following 
observation  about  Figure  4.  In  the  ideaHzed  Brownian  setting,  when  the  scaled  workload 
process  W  is  on  the  lower  ray  of  the  cone  boundary  and  Wi{t)   <  xi,  then  there  are 
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zero  scaled  customers  at  station  2  and  yet  station  2  is  not  idle.  This  apparent  paradox 
is  due  to  the  rescaling  that  occurs  when  passing  to  the  heavy  traffic  hniit.  In  the  actual 
queueing  system,  there  are  enough  customers  at  the  particular  station  to  avoid  idleness, 
but  when  looked  at  in  the  scaled  space  of  the  heavy  traffic  limit,  these  customers  vjinish. 
This  difference  may  prevent  the  workload  regulating  input  policy  to  achieve  the  desired 
throughput. 

However,  the  release  rule  can  be  adapted  to  the  actual  queueing  system  by  enlarging 
the  release  region.  There  are  two  ways  this  can  be  achieved.  The  first  way  is  to  slightly 
enlarge  the  region  on  the  inside  of  the  cone  boundary.  This  is  done  by  translating  the 
vertex  of  the  cone  from  (0, ...,  0)  to  (e, ...,  e).  This  translation  changes  the  hyperplcine  (8.2) 
to 

Y,aiW,it)  =  e(^a,).  (8.5) 

i-l  1=1 

This  treinslation,  which  may  be  negligible  in  scaled  space,  prevents  the  process  W  from 
straying  very  far  from  the  original  truncated  cone  bovmdary.  The  second  way  to  enlarge 
the  releeise  region  is  to  inflate  the  bounded  region  derived  in  Section  6  by  a  constant 
factor  K  >  1,  while  still  maintaining  the  relative  shape  of  the  bounded  region.  As  e  and 
K  increase,  the  servers  will  incur  less  idleness  but  the  queue  lengths  may  grow  as  a  result. 
The  workload  regulating  input  policy  sets  the  parameter  e  and  k  so  that  the  desired  output 
rate  A  is  achieved.  In  an  actual  queueing  system,  the  setting  of  e  and  k  will  depend  upon  a 
variety  of  factors,  including  the  amount  of  variability  in  the  queueing  system,  the  amiount 
of  time  customers  spend  at  non-bottleneck  stations,  the  network  topology,  and  the  extent 
to  which  the  network  is  heavily  loaded. 

Notice  that  the  input  policy  defined  in  this  section  has  been  described  in  terms  of 
the  seeded  workload  process  W.  Thus,  before  implementing  this  policy  in  an  actual 
queueing  system,  the  release  region  needs  to  be  expressed  in  terms  of  the  actual  (un- 
sealed) workload  process,  which  is  denoted  by  w{t)  =  {{wi{t),  ...,wj{t)),t  >  0}.  Then 
"^«(0  =  Jlk=i  ^^ikQki^)  for  i  =  1,...,/  and  t  >  0,  where  Qk{t)  is  the  actual  number  of 
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class  k  customers  in  the  system  at  time  t.  By  equation  (2.2)  it  follows  that 

XV  [Tit] 

Wi{t)  =      'y'  for  1  =  1, ...,  /    and  <  >  0.  (8.6) 

Since  the  scheduling  problem  was  solved  using  the  long-run  average  criterion,  the  time 
scaling  can  be  ignored,  and  replacing  Wi{t)  by  Wi{t)/y/n  in  the  inequedities  defining  the 
release  region  (8.4)  will  lead  to  an  implementable  release  policy  for  the  original  queueing 
network. 


9.  The  Workload  Balancing  Input  Heuristic 

The  customer  release  policy  described  in  the  last  section  allows  for  the  controller  to 
decide  when  to  release  the  next  customer  into  the  system,  but  not  to  choose  the  class  of 
the  entering  customer.  It  was  assumed  that  the  class  designations  of  entering  customers 
are  exogenously  chosen  so  that  q^  is  the  long-run  proportion  of  class  k  customers  released 
into  the  network.  In  this  section,  we  develop  a  workload  balancing  input  heuristic,  which 
is  based  on  insight  gained  from  the  solution  to  the  constrained  singular  control  problem, 
that  decides  which  class  of  customer  to  release  into  the  system.  This  scheme  appears  to 
improve  the  performance  of  the  scheduling  policy  and  is  guaranteed  to  keep  the  actual  mix 
of  released  customers  sufficiently  close  to  the  desired  mix  q. 

The  key  idea  behind  the  heuristic  is  the  observation  that  server  idleness  is  incurred  in 
the  idealized  Brownian  network  only  when  the  (/  —  1)— dimensional  workload  imbalance 
process  W  reaches  the  reflecting  boundary  derived  in  Section  6.  The  workload  imbalance 
process  stays  within  a  region  containing  the  origin,  but  when  it  reaches  the  reflecting 
boundary,  the  workload  becomes  too  imbalanced,  the  control  U  is  exerted,  and  at  least 
one  server  incurs  idleness. 

Thus,  server  idleness  would  be  reduced,  and  hence  system  performance  would  be  im- 
proved, if  the  workload  imbalance  process  could  be  discouraged  from  reaching  the  reflecting 
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boundaries.  Recall  by  equations  (3.13)  and  (7.2)  that 

K 
W,{t)  =  Y^  M,kZk{t)   for  z  =  1, ...,/-  1  and  <>  0.  (9.1) 

k=\ 

This  equation  relates  the  workload  imbaJajice  process  W  to  the  vector  queue  length  process 

Z.  Our  heuristic  will  release  the  customer  class  k  that  attempts  to  balsuice  the  workload 

imbalance  process  and  hence  avoid  server  idleness. 

Consider  the  actual  (unsealed)  queueing  system,  where  w,{t)  is  the  actual  workload 

imbalance  process  at  time  t  and  Qk(i)  is  the  actual  number  of  class  k  customers  in  the 

system  at  time  t;  then  Wi{t)  =  ^^=1  ^ikQk(t)  for  i  =  l,---,-^  —  1  and  t  >  0.    Suppose 

the  input  poHcy  derived  in  the  last  section  dictates  that  a  customer  is  to  be  released  into 

the  system  at  time  t.  There  are  two  steps  in  the  workload  balancing  input  heiiristic.  The 

first  step  checks  to  see  if  the  actual  mix  of  customers  already  released  into  the  network 

is  sufficiently  close  to  the  derived  mix  q.    Let  Nk{t)  denote  the  total  number  of  class  k 

customers  released  into  the  system  in  the  time  interval  [0,  t].  Let  E  —  {k  =  I,  ...,K\qk  >  0} 

be  the  set  of  possible  entering  classes;  these  classes  correspond  to  the  first  stage  of  some 

customer  type's  route.  Consider  the  constraints 

9jNj{t)  -  qkNk{t)  <  N*   for  all  ;,  k  e  E,  (9.2) 

where  N*  is  an  exogenously  specified  parameter  that  specifies  how  close  the  actual  entering 

class  mix  must  stay  to  the  target  mix  q. 

K  any  of  the  constraints  in  (9.2)  are  violated,  then  the  heuristic  releases  a  class  / 

customer,  where 

max  {qjNj{t)  -  qkNk{t)]  =  q,nNm{t)  -  qiN,{t).  (9.3) 

j,keE 

That  is,  we  release  the  customer  class  that  is  lagging  behind  the  farthest  from  its  desired 
target.  If  constraints  (9.2)  are  all  satisified,  then  the  actual  mix  of  released  customers  is 
suflRciently  close  to  the  desired  mix,  and  we  can  proceed  to  the  second  step  of  the  heuristic, 
which  attempts  to  balance  the  workload. 
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Let  Wi(t)  equal  the  time  average  value  of  Wi  over  the  time  interval  [0,t].  This  value 
can  be  easily  calculated  in  a  computer  simulation  model  of  a  queueing  network  or  in  an 
actual  factory  that  is  under  computer  control;  if  the  factory  is  not  under  computer  control, 
then  the  value  of  Wi(t)  may  be  recorded  periodically,  for  example  once  per  shift,  and  Wi{t) 
can  be  updated  accordingly.  Let  Wikit)  be  defined  by 

Wik{t)  =  w^{t)  +  Mik  for  z  =  1, ...,  J  -  1  and  k  ^  E.  (9.4) 

Thus,  Wiki^)  is  the  zth  component  of  the  workload  imbalance  process  that  would  result 
if  a  class  k  customer  were  released  into  the  system  at  time  t.  Step  two  of  the  workload 
balancing  input  heuristic  releases  a  class  /  customer  into  the  network,  where 

1=1  1=1 

Thus,  step  two  attempts  to  push  the  workload  imbalance  process  toward  its  long-run 

average  value,  which  presumably  will  not  be  close  to  the  reflecting  boundary.  Notice 
that  the  time-average  value  of  w  is  chosen  as  the  desired  target,  as  opposed  to  choosing 
the  origin  (that  is,  the  point  (0,  ...,0))  as  the  target.  This  is  because,  depending  on  the 
workload  profile  matrix  M  and  the  topology  of  the  network,  it  is  possible  that  the  origin 
will  not  be  a  particularly  desirable  value  of  u',  in  terms  of  avoiding  server  idleness. 

The  workload  balancing  input  heuristic  is  equally  applicable  to  multiclass  closed 
queueing  networks  (see  Section  4  of  Harrison  and  Wein  [8])  because  the  same  relationship 
holds  between  server  idleness  and  the  workload  imbalance  process.  In  a  closed  network,  a 
new  customer  is  released  into  the  network  whenever  a  customer  exits,  and  this  heuristic 
can  be  used  to  decide  which  class  of  customer  to  release.  This  heuristic  was  tested  on  a 
simulation  model  of  the  two-station  closed  network  example  in  Section  6  of  Harrison  and 
Wein  [8].  The  example  there  had  two  customer  types,  denoted  by  A  and  B,  and  the  desired 
input  mix  was  50-50.  The  workload  balancing  input  heiu-istic  offered  a  7.8%  improvement 
in  average  cycle  time  (from  54.9  to  50.6;  see  Table  I  of  [8])  over  deterministic  input  (releas- 
ing customers  in  the  order  ABABAB...),  while  maintaining  the  same  average  throughput 
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rate.  Similarly,  the  heuristic  offered  a  10.1%  reduction  in  average  cycle  time  (from  38.6  to 
34.7)  for  the  controllable  input  case  (see  Table  I  of  Wein  [21]).  For  both  of  these  cases, 
the  exogenous  parameter  TV*  was  chosen  to  guarantee  that  the  resulting  class  mix  wo\ild 
be  within  one-half  of  1%  of  the  desired  50-50  mix. 

Furthermore,  this  simple  heuristic  is  probably  effective  for  any  factory,  regardless  of 
the  timing  of  its  input  policy  (exogenous,  closed,  or  controllable)  or  priority  sequencing 
policy.  The  (/  —  1)— dimensional  workload  imbedance  process  offers  a  concise  and  effective 
measure  of  the  balance  of  work  tliroughout  a  complex  network,  and  its  crucial  relationship 
to  the  server  idleness  process  is  exploited  by  this  heuristic. 

10.  An  Example 

The  procedure  described  in  this  paper  will  be  illustrated  by  means  of  a  three-station 
example.  The  example  has  three  customer  types,  denoted  by  A,  B,  and  C,  and  the  specified 
product  mix  is  to  have  equal  numbers  of  aJl  three  types.  Table  I  describes  the  deterministic 
route  of  each  customer  type,  and  gives  the  mean  processing  time  (in  arbitrary  time  units) 
for  each  of  the  various  stages  of  service.  All  service  time  distributions  are  assumed  to  be 
exponential,  although  our  results  hold  for  any  service  time  distribution  with  finite  mean 
and  variance. 

MEAN 
CUSTOMER  SERVICE 

TYPE  ROUTE  TIMES 

A  3^  1  ^2  G.O   4.0    1.0 

B  1-^2-+ 3^  1^2  8.0   6.0    1.0   2.0    7.0 

C  2-^3^  1  ^3  4.0   9.0   4.0   2.0 

Since  each  customer  class  corresponds  to  a  type-stage  pair,  the  twelve  customer  classes 
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are  designated  (and  ordered  from  k  =  1,...,12)  by  (A1,A2,A3,B1,...,B5,C1,...,C4).  The 
12  X  12  routing  matrix  P  has  non-zero  entries  P12  =  P23  =  -^45  =  -Pse  —  Pei  —  Pjs  = 
Pg  jQ  =  Pjo  11  =  Pii,i2  =  1-  Calculation  of  the  3  x  12  workload  profile  matrix  M  yields 

/4     40     10      22204       44     0\ 
M=l     11     13     13     777     4      0      00,  (10.1) 

\6     0     0      1       1      1     0     0     11     11     2     2/ 

where  M,jt  is  the  expected  remaining  processing  time  at  station  i  for  a  class  k  customer 
until  that  customer  exits  the  network.  Since  9  =  (i00|0000|00  0)'^,  equation 
(2.12)  yields  Vi  =  V2  =  U3  =  6,  thus  implying  perfect  system  balance  by  equation  (2.13). 
Therefore,  the  p  values  can  be  factored  out  of  equation  (7.2),  as  mentioned  in  Section  4, 
and  the  2  x  12  workload  imbalance  profile  matrix  M  can  be  given  by 

.^.       /-2     4     0      9       1      1     2     0     -7      -7       2       -2\  . 

^^=[-5     1     1     12     12     6     7     7     -7     -11     -2     -2  j  "  ^^^'^^ 

The  exogenous  output  rate  is  .15  customers  per  tmit  of  time,  so  that,  by  (2.13),  pi  =  P2  = 
P3  =  .9,  and  the  system  parameter  value  of  n  =  100  can  be  chosen.  The  holding  costs  are 
cjt  =  1  for  A-  =  1, ...,  12,  so  that  the  objective  is  to  minimize  the  long-run  expected  average 
number  of  customers  in  the  system  (or  the  long-run  expected  average  cycle  time)  subject 
to  meeting  the  long-run  expected  average  output  rate  of  .15  customers  per  unit  time. 
The  dual  LP  (3.9)-(3.10)  can  be  expressed  as 

max     Wiit)ni{t)  +  W2{t)TT2it)  (10.3) 

subject  to    M'^irit)  <    e,  (10.4) 

where  M  is  given  in  (10.2).  This  problem  can  be  solved  graphically  for  all  values  of 
the  workload  imbalance  process  W,  and  the  solution  is  given  in  Table  II.  There  are  six 
extreme  points  of  the  constraint  set  (10.4),  and  thus  the  scaled  workload  process  W  lives  on 
a  polyhedral  cone  with  six  faces.  The  six  extreme  points  lead  to  the  workload  imbalance 
space  P^  being  partitioned  into  six  regions,  which  were  referred  to  in  Section  8  as  P/, 
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and  which  axe  niimbered  in  Table  II  for  future  reference.  For  each  of  the  six  regions, 
dynamic  reduced  costs  are  calculated  according  to  (3.11)  ajid  a  priority  sequencing  policy 
is  developed  according  to  rules  (7.1)-(7.3).  The  resulting  sequencing  policy  is  given  in 
Table  III. 

In  order  to  find  the  workload  regulating  input  policy,  a  numerical  solution  is  needed 
to  problem  (4.24)-(4.26).  The  objective  function  cost  h{W{t))  defined  in  (4.2)  is  found 
from  n*{t)  in  Table  II.  By  (2.14),  the  righthand  side  values  in  (4.25)  are  7i  =  72  =  73  =  1- 
By  (4.4),  the  drift  of  the  two-dimensional  Brownian  motion  process  B  is  (0,0)  and  the 
calculations  in  (2.7),  (2.17),  (2.18),  and  (4.3)  yield  the  covariance  matrix 

12.333      6.778 


a  = 


6.778      12.444 


(10.5) 


REGION  NUMBER 
AND  DESCRIPTION 


1:  t^i(0>0,W'2(0>0,|<{j;^<4; 
2:  W,it)>0,W,it)>0,^^<p^^<l; 
3:  Wi{t)  <  0,W2{t)  >  0; 
W,{t)<0,W2{t)^0; 

1^1(0  <0,H'2(0<0,^>1; 
4:  W,ii)<0,W2{t)<0,^<'Mll<i- 

5:  Wi{t)  =  0,W2it)  <0; 


W^i(0>0,H^2(^)<0,|^>-l; 
6:  Wi{t)>0,W2{t)  =  0\ 


►2(0 
i'2(0 


H^i(0>0,U2(0>0,^>4; 
Wyit)>0,W2ii)<0,p^<-l; 


DUAL 

SOLUTION 

7r*(0  =  -^,7r*(0  =  A 


7rJ(^)  =  -l,7r»(O  =  0 


^m  =  \^^2ii)  =  -l 


'^KO  =  ^>'^2(0  =  -^ 


TABLE  II.  Dual  Solution  as  a  Function  of  Workload  Imbalance 
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Using  a  mesh  size  oi  h  =  0.5  in  the  finite  difference  approximation,  the  solution  to  the 
LP  (6.1)-(6.4),  (6.8)  gives  the  reflecting  boundary  shown  in  Figure  5.  Superimposed  on  top 
of  this  boundary  is  the  partition  of  the  six  regions  i?y,  /  =  1, ...,  6,  corresponding  to  the  six 
dual  LP  solutions  described  in  Table  IL  The  workload  imbalance  process  is  uncontrolled 
when  it  is  strictly  within  the  reflecting  boundary,  and  when  it  reaches  the  boundary,  the 
process  is  pushed  back  in  (in  the  direction  of  the  arrows  in  Figure  5)  by  at  Iceist  one  of  the 
three  controls  (t/*,  C/2  >  ^3*)-  Notice  that  there  axe  places  on  the  boundary  where  more  than 
one  control  is  used  at  a  given  time.  In  particular,  U*  and  U^  are  used  at  states  (0.5,9.5), 
(1,9.5),  and  (1.5,10),  and  l/j  ^^^  ^3  ^^^  both  used  when  the  process  reaches  (3,0),  (3,0.5), 
(3,1),  and  (3,1.5). 


REGION 

STATION  1 

STATION  2 

STATION  3 

1 

B4  C3  A2  Bl 

A3  B5  B2  CI 

C4  B3  C2  Al 

2 

C3  A2  B4  Bl 

A3  CI  B5  B2 

C4  Al  C2  B3 

3 

A2  C3  B4  Bl 

A3  B5  CI  B2 

C4  Al  B3  C2 

4 

A2  C3  B4  Bl 

A3  B5  B2  CI 

C4  B3  Al  C2 

5 

B4  Bl  A2  C3 

A3  B5  B2  CI 

C4  B3  Al  C2 

6 

B4  Bl  C3  A2 

A3  B5  B2  CI 

C4  B3  Al  C2 

TABLE  III.  Priority  Sequencing  Policy 

As  mentioned  in  Section  8,  the  numerical  solution  to  the  constrained  singular  control 
problem  needs  to  be  transformed  into  an  explicit  expression  for  the  reflecting  boundary.  A 
piecewise  linear  boundary  consisting  of  13  segments  (see  Figure  6)  was  used  to  approximate 
the  reflecting  boundary  from  Figure  5.  The  release  region  (8.4)  for  this  example  is  given 
in  terms  of  the  actual  (unsealed)  workload  process  w  in  the  Appendix;  the  bounded  region 
Rb  is  expressed  in  the  Appendix  in  terms  of  the  workload  process  rather  than  the  workload 
imbalance  process.  Also,  recall  that  the  multipHcative  factor  k  appearing  in  the  Appendix 
inflates  the  bounded  region  in  Figure  6. 
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W2(t) 


-6-5-4-3-2-1  0         12  3  4  5 


FIGURE  5.  The  Reflecting  Boundary  in  the  Constrained  Singular  Control  Problem. 
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FIGURE  6.  A  Piecewise  Linear  Approximation  to  the  Reflecting  Boundary. 
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Using  the  example  network,  a  simulation  study  was  undertciken  to  compEire  the  perfor- 
mance of  the  scheduling  policies  proposed  in  Sections  8,  9,  and  10  agaiinst  several  conven- 
tional customer  release  and  priority  sequencing  policies.  Two  conventional  input  policies 
were  tested:  deterministic,  where  the  interarrival  times  axe  ail  constant,  and  closed  loop 
input,  where  the  total  number  of  customers  in  the  network  is  held  constant  at  A'^;  the 
latter  pohcy  is  abbreviated  by  CLOSED(N)  in  Table  IV.  For  both  of  these  input  rules, 
customers  were  released  into  the  network  in  the  order  ABCABCABC...  Both  input  pohcies 
were  tested  in  conjunction  with  two  priority  sequencing  rules:  first-in  first-out  (FIFO)  and 
the  shortest  expected  remaining  processing  time  rule  (SERPT),  where  priority  is  given  to 
the  customer  class  k  with  the  smallest  value  of  X^,_2  M,k. 

The  scheduling  policy  proposed  here  is  abbreviated  in  Table  IV  by  WR(e,  k), 
WBAL(A'^*),  and  DRC  (for  dynamic  reduced  costs),  where  the  workload  regulating  in- 
put policy  with  parameters  e  and  k  dictates  via  (8.4)  when  a  customer  is  to  be  released, 
the  workload  balancing  input  heuristic  with  pju-ameter  A'^*  dictates  the  class  of  customer 
to  be  released,  and  the  priority  sequencing  policy  is  defined  by  the  dynamic  indices  in  (7.1) 
and  (7.3). 


RULE 

RULE 

DETERMINISTIC 

FIFO 

DETERMINISTIC 

SERPT 

CL0SED(18) 

FIFO 

CLOSED(25) 

SERPT 

WR(1.7,1.5),WBAL(f )      DRC 


THROUGHPUT 
RATE 
(95%  C.Ll 
.149(±.0001) 
.149(±.0002) 
.149(±.0010) 
.149(±.0008) 
.149(±.000S) 


CYCLE 

TIME 
(95%  C.I.) 
144(±10.4) 
1S2(±15.7) 
120(±0.8) 
166(±1.1) 
S5.4(±l.l) 


TABLE  IV.  Comparison  of  Cycle  Times 

The  results  of  the  simulation  study  are  summarized  in  Table  IV.  Each  row  gives 
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statistics  for  a  particular  scheduling  policy,  which  is  a  combination  of  a  customer  release 
policy  and  a  priority  sequencing  rule.  For  each  policy  tested,  20  independent  runs  were 
made,  each  consisting  of  5000  customer  completions.  The  first  200  time  units  of  each  rim 
was  truncated  to  reduce  the  initial  bias.  The  third  and  fourth  coliunns  give  the  average 
throughput  rate  and  average  cycle  time,  respectively,  over  the  20  runs,  cdong  with  95% 
confidence  intervals.  The  parameters  N,  e,  and  k  were  chosen  so  that  all  scheduling  policies 
achieved  the  average  throughput  rate  of  .149  customer  completions  per  unit  time,  which 
corresponds  to  a  server  utilization  of  89.4%.  Recall  that  the  objective  is  to  minimize  the 
average  cycle  time  subject  to  a  given  average  throughput  rate. 

Referring  to  Table  IV,  it  is  seen  that  the  scheduling  policy  proposed  in  this  paper  easily 
outperforms  all  of  the  conventional  scheduling  rules.  The  policy  offers  a  28.8%  reduction 
in  average  cycle  time  over  the  (CLOSED, FIFO)  case,  which  was  its  closest  competitor. 
It  also  achieved  a  40.7%  reduction  relative  to  the  (DETERMINISTIC,FIFO)  case,  even 
though  it  is  well  known  (see  Whitt  [23],  for  example)  that  reducing  the  variability  in  the 
interarrival  times  of  a  queueing  network  will  lead  to  improved  performance. 

The  parameter  value  e  =  1.5  corresponds  to  15.0  units  of  unsealed  work,  which  in 
turn  roughly  corresponds  to  the  amount  of  work  for  each  server  that  is  embodied  in  two 
and  one-half  customers,  since  ii,  =  6  for  i  =  1,2,3.  The  parameter  value  «;  =  1.7  means 
that  the  bounded  region  in  Figure  6  was  enlarged  by  70%.  The  parameter  value  A'^*  =  ^ 
guaranteed  that  the  actual  entering  fraction  of  each  customer  class  was  within  .003  of  the 
specified  target  of  g^,.  =  .333;  in  the  simulation  study,  the  resulting  class  mix  was  virtually 
equal  to  the  target  mix. 

11.  Conclusions 

In  this  paper  we  have  considered  the  problem  of  how  to  dynamically  release  jobs 
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(when  and  which  class)  and  prioritize  jobs  in  a  multistation  multiclass  queueing  network 
with  general  service  time  distributions  and  a  general  routing  structure.  The  objective  weis 
to  minimize  the  long-run  expected  average  hnear  holding  costs  of  customers,  subject  to 
a  specified  average  input  mix  and  a  constraint  on  the  long-run  expected  average  output 
rate.  Under  baJemced  heavy  loading  conditions,  the  scheduling  problem  was  approximated 
by  a  Brownian  control  problem,  and  a  numerical  solution  to  the  workload  formulation  of 
the  control  problem  was  obtained. 

This  solution  was  then  interpreted  in  terms  of  the  original  queueing  system  in  order  to 
develop  an  effective  three- part  scheduling  policy.  The  first  part  is  the  workload  regulating 
input  policy,  which  releases  a  job  whenever  the  workload  process  enters  a  particular  region. 
The  second  part  is  the  workload  balancing  input  heuristic,  which  releeises  the  customer 
class  that  will  best  balance  the  workload  among  the  various  bottleneck  stations.  The  third 
part  is  the  priority  sequencing  policy,  which  assigns  dynamic  indices  (based  on  dynamic 
reduced  costs  from  a  linear  program)  to  each  customer  class.  A  computational  study  was 
performed  that  exhibited  the  policy's  effectiveness. 

Two  related  research  topics  are  to  prove  a  weak  convergence  result  for  the  finite 
difference  approximation  used  in  Section  5,  and  to  develop  an  efficient  algorithm  to  solve 
the  large  linear  program  posed  in  Section  6.  Also,  the  close  relationship  between  the 
problem  addressed  in  this  paper  and  the  problem  of  priority  sequencing  in  a  multistation 
multiclass  closed  queueing  network  requires  further  investigation.  Finally,  more  numerical 
studies  need  to  be  performed  to  better  understand  the  behavior  and  robustness  of  the 
proposed  scheduling  policy. 

Appendix 

Let  the  regions  il/,  /  =  1, ...,  6  be  defined  as  in  the  first  column  of  Table  II  (with  W{t) 
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replaced  by  u'(t)).    The  release  region  (8.4)  for  the  example  in  Section  10  is  to  release  a 
customer  whenever  the  workload  process  w(t)  enters  ni</<6i?/,  where  Rj  =  RfDSf,  and: 

Si  :    wi{t)  -  4u'2(0  +  42ti'3(0  <  390e,  and 

Wi{t)-W3{t)  <  30k. 

^2  :    -W2(t)  +  13w3it)  <  120e, 
wi(t)  —  W3(t)  <  35k, 
2wi{t)  +  W2{t)  -  3w3{t)  <  140k,  and 

-Wi{t)  +  2u>2(<)  -  lV3{t)  <  ISOk. 

53  :    I39wi(t)  -  18w2{t)  -  44t('3(0  <    770e, 

W2{t)  —  W3{t)  <  95k,  and 
-wi{t)  +  W3{t)  <  0; 

or 
lZ9wi{t)  -  I8w2{t)  -  44u'3(0  <    770e,  and 
-Uwi{t)  +  5w2{t)  +  6w3{t)  <  300k. 

54  :    ll«'i(0-4u)3(0  <    70e, 

— u'i(i)  +  u'3(i)  <  55k,  and 
wi{t)  -W2{t)  <  25k. 

55  :    W2{t)  <    lOe,  and 

Wi{t)  -W2{t)  <  15k. 

56  :    -wi{t)  +  4u'2(0  +  2u'3(f )  <  50e,  and 

Wi(t)  —  W2{t)  <  15k; 
or 

-wiit)  +  4w2{t)  +  2iV3{f)  <50e, 
-W2{t)  +  iC3{t)  <  0,  and 
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Wi{t)-W3{t)  <  30/c. 
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